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Abstract 


An  algorithm  is  presented  to  simulate  fluid  dynamics  on  a  three  qubit  type  II 
quantum  computer:  a  lattice  of  small  quantum  computers  that  communicate  classical 
information.  The  algorithm  presented  is  called  a  three  qubit  factorized  quantum  lattice 
gas  algorithm.  It  is  modeled  after  classical  lattice  gas  algorithms  which  move  virtual 
particles  along  an  imaginary  lattice  and  change  the  particles’  momentums  using  collision 
rules  when  they  meet  at  a  lattice  node.  Instead  of  moving  particles,  the  quantum 
algorithm  presented  here  moves  probabilities,  which  interact  via  a  unitary  collision 
operator.  Probabilities  are  determined  using  ensemble  measurement  and  are  moved  with 
classical  communications  channels.  The  lattice  node  spacing  is  defined  to  be  a 
microscopic  scale  length.  A  mesoscopic  governing  equation  for  the  lattice  is  derived  for 
the  most  general  three  qubit  collision  operator  which  preserves  particle  number.  In  the 
continuum  limit  of  the  lattice,  a  governing  macroscopic  partial  differential  equation — the 
diffusion  equation — is  derived  for  a  particular  collision  operator  using  a  Chapman- 
Enskog  expansion.  A  numerical  simulation  of  the  algorithm  is  carried  out  on  a 
conventional  desktop  computer  and  compared  to  the  analytic  solution  of  the  diffusion 
equation.  The  simulation  agrees  very  well  with  the  known  solution. 
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TYPE  II  QUANTUM  COMPUTING  ALGORITHM 
FOR  COMPUTATIONAL  FLUID  DYNAMICS 


1  Introduction 

1.1  Overview 

In  1982  Richard  Feynman  proposed  building  a  computer  based  on  quantum 
mechanical  principles  to  efficiently  simulate  quantum  systems.  In  the  two  decades  since, 
significant  progress  has  been  made  both  theoretically  and  experimentally  towards  this 
end.  Following  Feynman’s  vision,  in  2002  the  Air  Force  Research  Laboratory  and  the 
Air  Force  Office  of  Scientific  Research  established  a  basic  research  theme  called 
Quantum  Computation  for  Physical  Modeling  [1],  The  goal  of  this  project  is  to  explore 
quantum  algorithms  and  practical  quantum  computers  to  model  dynamic  physical 
systems  with  an  exponential  increase  in  computational  efficiency. 

This  thesis  supports  this  goal  by  extending  an  algorithm  designed  to  model  fluid 
dynamics  using  a  lattice  of  interacting  quantum  systems.  This  algorithm  was  used  by 
Yepez  to  investigate  Navier-Stokes  equations  of  fluid  dynamics  [2-5],  the  diffusion 
equation  [6-9],  and  the  Burgers  equation  [10-12],  It  is  called  the  Factorized  Quantum 
Lattice  Gas  Algorithm  (FQLGA),  and  it  derives  its  name  from  algorithms  written  to 
model  fluid  dynamics  on  classical  computers.  These  classical  algorithms  model  particles 
on  a  spatial  lattice.  The  particles  move  from  lattice  node  to  lattice  node  in  discrete  time 
steps.  They  include  classical  collision  rules  to  control  how  particles  interact  when  they 
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meet  at  a  lattice  node,  where  the  distance  between  nodes  is  defined  to  be  a  microscopic 
scale  length.  Since  the  algorithms  consists  of  a  gas  of  virtual  particles  moving  on  a 
discrete  lattice,  it  is  called  a  Lattice  Gas  Algorithm  (LGA)  by  the  fluid  dynamics 
community,  or  a  Classical  Lattice  Gas  Algorithm  (CLGA)  by  the  quantum  computing 
community  to  distinguish  it  from  similar  quantum  lattice  gas  algorithms.  One  may  use  an 
ensemble  of  (quantum  or  classical)  lattice  gases  to  develop  a  finite  difference  equation 
known  as  the  mesoscopic  lattice  Boltzmann  equation.  From  this  equation  it  is  possible  to 
develop  a  macroscopic  effective  field  theory  in  the  continuum  limit  of  the  lattice, 
essentially  by  taking  the  Taylor  series  expansion  of  the  lattice  Boltzmann  equation 
around  local  equilibrium  in  what  is  known  as  a  Chapman-Enskog  expansion. 

The  FQLGA  is  said  to  be  “factorized”  because  it  is  designed  to  run  on  a  quantum 
computer  that  is  not  fully  coherent — that  is,  on  a  computer  that  is  factorized  into  many 
smaller  quantum  computers  communicating  with  classical  information:  a  type  II  quantum 
computer.  This  sort  of  computer  is  interesting  because  a  prototype  already  exists  which 
has  run  fluid  dynamics  simulations  [6,  7,  11],  The  algorithm  developed  by  Yepez  [2-12] 
uses  a  lattice  of  two  quantum  bit  (qubit)  computers  to  perform  the  computations.  Instead 
of  moving  particles  across  a  one  dimensional  lattice  probabilities  are  moved,  and  instead 
of  using  collision  rules  to  govern  interactions  a  unitary  operator  called  the  collision 
operator  is  used.  This  thesis  extends  this  algorithm  to  run  on  a  three  qubit  type  II 
quantum  computer.  The  main  contributions  presented  in  this  paper  are  the  computation 
of  the  lattice  Boltzmann  equations  for  the  most  general  three  qubit  collision  operator  that 
conserves  particle  number,  and  the  derivation  of  the  diffusion  equation  as  an  effective 
field  theory  for  a  more  specific  collision  operator.  In  addition,  numerical  simulations 
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comparing  the  diffusion  equation  FQLGA  simulation  and  the  analytic  solution  of  this 
partial  differential  equation  are  presented. 

This  thesis  is  meant  to  be  accessible  to  someone  with  a  reasonable  background  in 
quantum  physics,  but  with  little  exposure  to  the  subjects  of  quantum  computing  or  fluid 
dynamics.  As  such,  I  will  introduce  some  of  the  basics  in  each  of  these  subjects  before 
discussing  the  FQLGA. 

1.2  Organization 

This  thesis  begins  with  a  summary  of  quantum  computing  in  Section  2.  This 
section  is  aimed  towards  those  who  have  studied  quantum  mechanics  but  have  had  little 
exposure  to  quantum  computing.  It  briefly  discusses  what  a  qubit  is,  how  quantum  logic 
gates  perform  computations,  and  how  qubit  measurement  affects  the  type  of  information 
one  may  obtain  from  it.  I  attempt  to  make  this  discussion  easier  to  follow  by  using 
analogies  with  the  more  familiar  classical  bits  and  classical  logic  gates. 

Following  this,  the  main  categories  of  quantum  algorithms  developed  thus  far  are 
reviewed  to  give  readers  an  idea  of  where  the  FQLGA  fits  in  the  world  of  quantum 
computing.  Subsequently,  several  types  of  quantum  computers  in  development  are 
discussed,  along  with  the  various  challenges  associated  with  constructing  each  kind  of 
computer.  The  focus  of  this  section  is  on  Nuclear  Magnetic  Resonance  (NMR)  quantum 
computers  because  the  most  advanced  quantum  computer  prototype  to  date  uses  NMR 
technology,  and  because  these  machines  are  best  suited  for  the  quantum  algorithm 
developed  in  this  paper. 
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Since  the  FQLGA  models  fluid  dynamics,  Section  3  introduces  the  basics  of 
Navier-Stokes  fluid  dynamics  as  well  as  classical  lattice  gas  and  lattice  Boltzmann 
algorithms  used  to  simulate  these  fluids.  These  classical  algorithms  can  be  used  as  an 
analogy  to  better  understand  their  quantum  counterparts.  A  number  of  important 
macroscopic  dimensionless  parameters  used  to  characterize  Navier-Stokes  fluids  are 
listed  and  drawn  on  during  an  explanation  of  diffusive  ordering.  The  lattice  Boltzmann 
equation,  which  comes  from  the  more  familiar  Boltzmann  equation,  is  also  introduced. 

Following  this  brief  introduction  to  quantum  computing  and  classical  lattice  gas 
algorithms,  Yepez’s  factorized  quantum  lattice  gas  algorithm  is  introduced  in  Section  4. 
The  details  of  this  algorithm  are  laid  out  and  the  unitary  collision  operator  is  introduced. 
Following  this  is  an  explanation  of  how  the  entire  algorithm  can  be  contained  in  a  single 
equation,  the  quantum  lattice  Boltzmann  equation  (QLBE).  Next,  the  local  equilibrium 
probabilities  are  derived  and  subsequently  used  in  the  Chapman-Enskog  expansion  of  the 
QLBE  to  derive  the  governing  effective  field  theory  of  the  lattice  in  the  continuum  limit. 
This  governing  equation  turns  out  to  be  the  one  dimensional  Burger’s  equation,  which  is 
a  second  order  nonlinear  partial  differential  equation  used  to  model  turbulence  and  shock 
formation  in  inelastic  gases. 

In  Section  5  the  three  qubit  FQLGA  I  have  developed  is  introduced  and  the  most 
general  three  qubit  collision  operator  that  conserves  particle  number  is  derived.  The 
quantum  lattice  Boltzmann  equation  is  obtained  using  this  operator. 

Section  6  introduces  a  more  specific  collision  operator  which  yields  the  diffusion 
equation  as  the  algorithm’s  continuum  limit  governing  partial  differential  equation.  A 
complete  derivation  of  this  equation  is  given  starting  from  the  QLBE  in  section  6.1.  In 
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section  6.2  the  results  of  a  numerical  simulation  carried  out  on  a  conventional  computer 
are  presented,  and  compared  to  the  analytic  solution  of  the  diffusion  equation.  Finally,  a 
discussion  of  the  error  and  of  the  convergence  properties  of  the  algorithm  as  compared  to 
the  analytic  solution  is  included. 

2  Quantum  Computing  Summary 

Quantum  computing  was  first  proposed  by  Richard  Feynman  in  1982  [13,  14]. 

He  noted  that  there  were  certain  difficulties  in  simulating  quantum  mechanical  systems 
on  classical  computers  due  to  the  exponential  growth  of  the  problem  with  system 
complexity.  He  suggested  developing  a  computer  based  on  the  principles  of  quantum 
mechanics  to  overcome  these  difficulties.  In  1985  David  Deutsch  expanded  on  this  idea 
while  trying  to  use  the  laws  of  physics  to  derive  a  stronger  version  of  the  Church-Turing 
thesis.  This  thesis  states  that  the  Turing  model  of  computation  is  at  least  as  efficient  as 
any  other  model  of  computation,  in  the  sense  that  if  one  computational  model  can  solve  a 
problem  in  time  polynomial  to  the  size  of  a  problem,  then  a  probabilistic  Turing  machine 
can  too  [15].  Since  the  laws  of  physics  are  ultimately  quantum  mechanical,  this  led 
Deutsch  to  develop  the  modem  concept  of  quantum  computers,  which  are  able  to 
efficiently  solve  problems  that  are  believed  to  have  no  efficient  solutions  on  classical 
computers  and  Turing  machines. 

Deutsch  developed  a  simple  algorithm  that  suggested  quantum  computers  would 
indeed  have  more  computational  power  than  a  classical  Turing  machine  (classical 
computer).  Over  the  next  decade  additional  algorithms  were  developed  culminating  in 
1994  with  Peter  Shor’s  factoring  and  discrete  logarithm  algorithms  [16],  and  in  1997  with 
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Lov  Grover’s  search  algorithm  [17].  However,  despite  the  immense  progress  made  in 
developing  these  algorithms,  physicists  have  so  far  managed  only  modest  advancements 
in  developing  physical  quantum  computers. 

Sections  2. 1  through  2.3  are  meant  to  give  a  brief  overview  on  how  quantum 
computing  works. 

2.1  Quantum  bits 

In  a  quantum  computer,  qubits  replace  classical  bits.  Qubits  are  analogous  to 
classical  bits  in  that  when  read  (measured)  they  can  only  be  0  or  1 .  However,  before 
measurement — when  a  computation  is  being  performed — a  qubit  can  be  in  a 
superposition  of  0  and  1  states.  Equation  (2.1)  shows  the  most  general  state  of  a  qubit 
where  1 1)  and  1 0)  are  the  basis  states  for  an  arbitrary  two  level  quantum  system  (for 
example  spin  up  and  down  in  a  spin  half  particle). 

|  q)  =  el *  sin(<9)  1 1)  +  cos(0)  1 0)  (2.1) 


2.1.1  Quantum  gates 

To  carry  out  a  quantum  algorithm,  quantum  computers  perform  unitary  operations 
on  qubits.  This  is  analogous  to  a  classical  algorithm  being  composed  of  ‘gates’  (AND, 
OR,  NOT,  XOR,  etc.)  that  act  on  classical  bits.  A  straightforward  example  of  this 
analogy  is  the  classical  exclusive  or  (XOR)  gate  and  the  quantum  controlled  not  (CNOT) 
gate.  The  XOR  gate  is  shown  in  Figure  la. 

The  values  ci  and  C2  represent  classical  bits  flowing  down  a  wire  (black  line)  from 
left  to  right.  They  encounter  the  XOR  gate  and  undergo  modulo  two  addition  (denoted 
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a 


b 


Figure  1.  a)  An  XOR  gate,  b)  A  modified  XOR  gate  that  preserves  the  bit  Ci. 


by  ®)  so  that  the  XOR  gate  output  is  0  ®  0  =  0,  0®  1  =  1,  1®0=  1,  and  1®  1  =0.  This 
is  described  in  the  truth  table  given  in  Table  la.  In  Figure  lb  we  have  created  a  modified 
gate  that  preserves  the  bit  ci.  Thus,  this  gate  has  an  equal  number  of  input  and  output 
bits.  All  quantum  gates  share  this  property  since  the  number  of  basis  vectors  used  to 
represent  |  y/)  remains  constant.  Two  electrons  will  always  be  spin  up  or  down  (or  both) 

no  matter  what  unitary  operations  one  performs  on  them,  so  it  is  impossible  to  create  a 
quantum  gate  with  fewer  outputs  than  inputs.  Therefore,  the  modified  XOR  gate  shown 
in  Figure  lb  will  be  a  better  analogy  to  the  CNOT  gate,  as  we  shall  now  see. 


Table  1.  a)  Truth  table  for  Figure  la.  b)  Truth  table  for  Figure  lb.  Inputs  are  left  of  the  gray  bar,  outputs 
are  right. 
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Figure  2  shows  a  quantum  CNOT  gate.  In  this  diagram,  straight  lines  represent 
qubit  states,  and  time  flows  from  left  to  right.  Again,  the  symbol  ®  denotes  modulo  two 
addition.  A  truth  table  for  this  gate  is  shown  in  Table  2.  Note  that  \q2)  passes 
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unchanged  as  long  as  \qt )  is  |0) ,  and  | q2)  is  changed  if  |g,)  is  |l) .  Thus  \qt)  is  called 
the  control  qubit  and  the  entire  gate  is  called  a  control  not  gate.  The  truth  table  given  in 
Table  2  is  the  same  as  Table  lb,  so  that  as  long  as  \qt)  and  \q2)  are  definitely  in  either 

state  1 1)  or  1 0)  the  gate  acts  as  a  modified  classical  XOR  gate.  However,  if  either  qubit 

is  in  a  superposition  of  1 1)  or  1 0)  states,  then  the  outputs  will  also  be  in  a  superposition 

of  states.  This  ability  of  quantum  bits  to  be  a  superposition  of  ones  and  zeros  at  the  same 
time  is  what  distinguishes  quantum  and  classical  computing. 


M 


Figure  2.  Controlled  not  gate. 

Table  2.  CNOT  gate  truth  table.  Inputs  are  left  of  the  gray  bar;  outputs  are  right. 
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The  state  of  the  entire  system  in  Figure  2  is  |  y/)  =  |  qx )  ®  |  q2 )  =  |  qx  q2) .  Note  that 
since  |</,)  acts  as  a  control,  it  is  possible  for  the  qubits  to  become  entangled.  For 
instance,  if  the  input  state  of  the  CNOT  gate  is  |  y/)  =  a  1 1  0)  +  b  1 0  0) ,  which  can  be 
factored  into  [a  1 1)  +  b 1 0))  0 1 0) ,  the  output  state  will  be  |  y/'j  =  a  1 1  1^  +  b 1 0  0^ ,  which 
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cannot  be  factored.  Before  application  of  the  CNOT  gate,  the  state  of  the  second  qubit 
was  independent  of  the  first  qubit,  allowing  |  i/z'j  to  be  factored.  However,  after 
application  of  CNOT,  the  state  of  qubit  two  was  directly  dependent  on  the  state  of  qubit 
one  and  |i//)  could  not  be  factored.  Thus,  after  the  CNOT  gate,  the  measurement  of  one 


qubit  will  immediately  determine  the  state  of  the  other  and  we  say  the  qubits  are 
entangled. 

The  CNOT  gate,  like  all  quantum  gates,  is  a  unitary  operation,  and  can  be 


described  in  matrix  form.  If  we  choose  the  following  to  be  our  basis 
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then  using  Table  2,  the  CNOT  matrix  is 


Ucnot 


^0  10  0^ 
10  0  0 
0  0  10 
v0  0  0  1, 


(2.2) 


(2.3) 


Returning  to  the  problem  discussed  earlier,  we  can  see  that  this  formalism  gives  the  same 
results. 


(a|l  0>+6|0  0»  =  (a|l  l)+6|0  0» 
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One  should  note  that  quantum  computers  have  an  exponential  increase  in 
computational  power  as  one  adds  more  qubits  [15].  This  is  because  each  additional  qubit 
doubles  the  number  of  Hilbert  space  dimensions  so  that  it  has  2s  dimensions,  where  B  is 
the  number  of  qubits.  For  instance,  if  we  have  three  qubits  there  are  23  dimensions  with 
one  choice  of  basis  being 
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For  purposes  of  comparison,  suppose  that  we  had  a  classical  computer  that  used 
three  parallel  lines,  each  of  which  simultaneously  transmitted  a  bit  to  a  central  processor. 
Then  at  any  given  time  the  processor  can  compute  using  three  bits,  which  we  will  define 
to  be  a  byte.  So  the  classical  computer  can  compute  with  one  byte  at  a  time.  In  contrast, 
a  similar  quantum  computer  with  a  three  qubit  memory — all  of  which  may  be  in  a 
superposition  of  1 1)  and  1 0)  states — can  compute  with  all  eight  bytes  simultaneously.  Of 
course  the  number  of  bytes  the  quantum  computer  can  handle  at  one  time  rises 
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exponentially  with  additional  bits.  Thus,  for  some  problems  the  quantum  computer  is 
exponentially  more  powerful  than  a  classical  computer. 


2.1.2  Measurement 

After  a  calculation  is  completed,  one  must  measure  the  quantum  bits  to  get  an 
answer.  Of  course,  if  the  output  is  in  a  superposition  of  states,  as  it  is  in  equation  (2.6) 
below,  then  one  will  get  a  random  answer  weighted  by  the  coefficients  in  front  of  each 
state.  For  the  equation  below,  the  computer  will  produce  the  binary  output  1,  1  with 

I  |2 

probability  \a\  . 


\y/)  =  a\\  1)  +  Z>|1  0)  +  c|0l)  +  </|00)  (2.6) 

To  avoid  the  embarrassment  of  getting  different  answers  each  time,  most  quantum 
algorithms  include  steps  to  make  the  coefficients  of  the  calculated  incorrect  answers  go  to 
zero.  However,  the  binary  ones  and  zeros  are  not  the  only  way  one  can  code  information; 
it  can  also  be  saved  in  the  magnitude  of  the  basis  coefficients.  To  get  this  information, 
one  must  either  perform  the  computation  many  times  on  the  same  computer  and  average 
the  measured  results,  or  perform  the  same  computation  on  many  identical  quantum 
computers  and  average  these  measured  results.  The  second  method  is  called  ensemble 


N  quantum  computers 

\y/)  =  a\\  1)  +  6|1  0)  +  c|0l)  +  <f|00) 
\y/)  =  a\\ l)  +  6|l  0)  +  c|0l)  +  J|00) 
\y/)  =  a\\ l)  +  6|l  0)  +  c|0l)  +  J|00) 
\y/)  =  a\\ l)  +  6|l  0)  +  c|0l)  +  J|00) 
\y/)  =  a\\ l)  +  6|l  0)  +  c|0l)  +  J|00) 


X 


j 


#  of  computers  in  1 1 1) 

N 

#  of  computers  in  |10) 

N 


etc... 


Figure  3.  Ensemble  measurement  averages  the  measurement  results  of  N  identical  quantum  computers  to 
obtain  the  magnitude  of  basis  coefficients.  The  symbol  with  the  arrow  in  the  figure  above  is  used  in 
quantum  computing  literature  to  signify  the  measurement  of  a  quantum  system. 
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measurement  and  is  how  information  is  extracted  in  a  Nuclear  Magnetic  Resonance 
(NMR)  quantum  computer.  The  FQLGA  takes  advantage  of  the  ability  of  NMR 
machines  to  do  ensemble  measurement,  so  these  computers  be  will  discussed  in  more 
detail  in  section  2.3. 

For  convenience,  one  may  represent  the  results  of  an  ensemble  measurement 
using  projectors  or  matrices.  For  instance,  suppose  we  are  interested  in  measuring  the 
probability  of  finding  the  second  qubit  in  Figure  3  in  the  state  1 1) .  Then  one  can  write 


p(M=\l))=\a\2 +\c\2 =|(11k)|2+|(01k>|2 =(H(I11)(11I+I01)(01l)k)-  (2-y) 

The  far  right  hand  side  of  (2.7)  indicates  that  the  probability  of  the  second  qubit  being  1 1) 
is  equal  to  (y  |  n2 1 y/) ,  where  n2  =  1 1  l)(l  1 1  + 1 0 1) (0 1 1  is  called  the  number  operator  and  is 
defined  to  be  the  sum  of  those  projectors  whose  second  qubit  is  |l) .  Using  the  basis 
given  in  (2.2),  we  can  rewrite  the  number  operator  in  matrix  notation: 


«2=|11)(11|  +  |01){01| 


O  0  0  0^ 
0  0  0  0 
0  0  10 
v0  0  0  oy 


(2.8) 


Similarly,  the  probability  of  finding  the  first  qubit  in  1 1)  is  equal  to  (y 1  h,  |  if/)  where 


«,=|11)(11|  +  |10)(10| 


^1  0  0  0^ 
0  10  0 
0  0  0  0 
v0  0  0  oy 


(2.9) 
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Note  that  h,  and  n2  are  not  unitary  and  therefore  do  not  represent  operations  that  can  be 

performed  by  a  quantum  computer.  Rather,  they  are  convenient  notations  allowing  one 
to  predict  the  results  of  an  ensemble  measurement. 


2.2  Quantum  algorithms 

A  quantum  algorithm  consists  of  a  series  of  unitary  transformations  performed  on 
qubits  followed  by  a  measurement  designed  to  perform  a  computation.  To  this  date,  the 
types  of  quantum  algorithms  developed  generally  fall  into  three  categories:  Fourier 
transform,  search,  and  simulation  algorithms  [15]. 

The  Fourier  transform  algorithm  is  the  backbone  of  Shor’s  factoring  and  discrete 
logarithm  algorithms,  and  it  involves  taking  the  Fourier  transform  of  a  set  of  numbers: 

{ x0,...  x  }  to  get  a  new  set:  { y(P...  y  „ }.  Suppose  one  prepares  a  state, 

2"-l 

W)=Hxj\j)  ,  so  that  the  coefficients  jc  .  of  the  basis  states  are  the  numbers  one  wishes 

j= 0 

to  transform.  Then  one  may  define  a  unitary  transformation  such  that 

i»->4r  ly**'2- 1*>  •  <2'io> 

V2  k= 0 

If  this  transformation  is  performed  on  |  y/) ,  we  see  that  the  new  coefficients  are  the 


Fourier  transformed  set  { y0,...  y  },  which  we  wanted. 


2-1 


7=0 


k= 0 


2n-\ 


1  y  e2*wrx 

V2" 


7=0 


2W-1 


k)=Hyk\k) 


(2.11) 


k=0 


Of  course  one  cannot  simply  read  off  the  coefficients  jy.  If  one  tried,  the  wave 
function  would  collapse  into  a  random  collection  of  bits.  It  takes  an  additional  amount  of 
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cleverness  and  a  few  more  quantum  logic  gates  to  get  useful  information  from  this 
transformation.  Nevertheless,  this  algorithm  can  complete  the  transform  in  about  n 2  steps 
as  opposed  to  the  classical  n2n  steps  (for  2"  numbers)  [16] — an  exponential  speedup! 

Grover’s  search  algorithm  is  an  example  of  the  second  kind  of  quantum 
algorithm — the  search  algorithm.  Grover’s  algorithm  is  designed  to  search  a  space  of 
size  n,  looking  for  an  element  in  it  with  some  desirable  attributes,  with  no  information 
about  the  structure  of  the  space.  Classically,  this  problem  requires  about  n  steps  while 

the  quantum  algorithm  can  accomplish  it  in  about  4n  steps  [17]. 

Finally,  simulation  algorithms  can  be  used  to  model  physical  (typically  quantum) 
systems.  Quantum  computers  are  ideal  for  this  task  because  their  Hilbert  space  increases 
exponentially  with  the  number  of  qubits  involved.  If  the  system  we  are  attempting  to 
model  is  quantum  mechanical  and  has  n  components,  then  in  general  it  takes  c"  bits  of 
memory  on  a  classical  computer  to  model  it,  where  c  is  some  constant  associated  with  the 
details  of  the  system.  On  the  other  hand,  a  quantum  computer  only  requires  k  n  qubits  to 
model  the  system,  where  again  A:  is  a  constant  that  depends  on  the  system  [15].  Though 
simulation  algorithms  intended  for  quantum  computers  are  typically  designed  to  model 
quantum  systems,  they  can  also  model  classical  systems.  The  Factorized  Quantum 
Lattice  Gas  Algorithm  is  an  example  of  a  quantum  algorithm  designed  to  model  a 
classical  system. 

2.3  Physical  quantum  computers 

Due  to  the  enormous  challenge  associated  with  isolating  and  precisely  controlling 
single  particles,  quantum  computers  are  currently  incapable  of  rivaling  their  classical 
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counterparts.  Quantum  computers  come  in  two  varieties  called  type  I  and  type  II.  Type  I 
machines  are  ‘pure’  quantum  computers  and  utilize  a  number  of  qubits,  each  of  which 
can  be  entangled  with  any  other  using  an  arbitrary  unitary  transformation.  Type  II 
machines  are  not  as  powerful  but  are  easier  to  create  in  practice.  They  consist  of  a 
number  of  small  type  I  quantum  computers  (called  nodes)  with  as  few  as  two  qubits  in 
each,  connected  by  classical  communications  channels  carrying  bits  instead  of  qubits 
[18].  Figure  4  shows  a  simple  diagram  of  a  type-II  computer. 


Claulcal 

Commniilcatloiu 

Chuual 


Figure  4.  Type-II  quantum  computer. 

The  four  most  developed  technologies  for  quantum  computing  are  optical 
techniques,  ion  traps,  neutral  atom  traps,  and  nuclear  magnetic  resonance.  One  of  the 
most  significant  problems  for  each  of  these  technologies  is  decoherence.  Decoherence  is 
the  uncontrolled  entanglement  of  a  system  with  its  environment,  destroying  the 
superposition  of  qubit  states  within  the  system  and  losing  the  information  it  contains  [19]. 
The  time  it  takes  for  this  process  to  occur,  called  the  decoherence  time,  is  very  short  in 
most  systems  and  therefore  limits  the  number  of  operations  that  can  be  performed  on  a 
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given  number  of  qubits.  Table  3  lists  the  various  decoherence  times  (  tq)  of  some  of  the 
systems  under  investigation,  along  with  the  time  it  takes  to  perform  an  operation  ( rop)  on 
the  system  [15].  This  gives  a  general  idea  of  the  number  of  operations  that  can  be 
performed  on  the  system  ( nop )  before  quantum  information  is  lost. 


Table  3.  Estimates  for  the  decoherence  time  Tq,  operation  time  rop,  and  maximum  number  of  operations  nop 
for  quantum  computer  candidates  [15]. 


System 

rQ  (sec) 

Top  (sec) 

N  op 

Nuclear  spin 

1  (T  -10* 

1 0"3  - 1 0’6 

10-  -1014 

Electron  spin 

IT3 

10-7 

104 

Ion  trap  (ln+) 

io- 

10-14 

1013 

Electron  -  Au 

10^ 

IQ-14 

10e 

Electron  -  GaAs 

1(T“ 

IQ-1 3 

103 

Quantum  dot 

10^ 

10’5 

103 

Optical  cavity 

10^ 

10-14 

10“ 

Microwave  cavity 

10': 

10-4 

104 

Given  the  long  decoherence  times  of  nuclear  spins,  it  is  not  surprising  that  the 
most  advanced  quantum  computers  rely  on  encoding  quantum  information  in  atoms  with 
spin  half  nuclei  using  NMR  technology.  This  is  done  by  placing  a  liquid  sample  in  a 
magnetic  field  around  10  T,  splitting  the  nuclear  spin  energy  levels  with  a  sort  of  Zeeman 
shift.  Radio  frequency  pulses  can  then  be  used  to  manipulate  the  nuclear  states. 


Figure  5.  The  molecule  used  by  Pravia  et  al  [6]  in  their  implementation  of  a  type  II  NMR  quantum 
computer  was  13C-Chloroform,  with  hydrogen  and  carbon  13  nuclear  spins  serving  as  the  qubits.  The 
energy  levels  of  the  nuclear  spin  states  are  split  using  a  strong  magnetic  field. 
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A  basic  description  of  a  type-II  NMR  quantum  computer  is  shown  in  Figure  5  and 
Figure  6.  The  nucleons  that  act  as  qubits  are  in  the  molecules  that  make  up  a  liquid 
sample — in  effect  each  molecule  is  a  small  quantum  computer.  Radio  frequency  pulses 
are  used  to  perform  unitary  transformations  on  the  qubits  in  the  sample.  Since  the  sample 
contains  some  1023  identical  molecules,  using  a  NMR  quantum  computer  amounts  to 
performing  the  same  calculation  on  10  quantum  computers.  Therefore,  when  one 
measures  the  state  of  a  particular  qubit,  one  does  it  for  the  entire  sample  and  gets  an 
average  value  of  the  state.  This  is  an  example  of  an  ensemble  measurement  discussed  in 
section  2.2. 


Gradient  Coil 


t 

Liquid  Sample  Separated 
into  Spatial  Nodes 


RF  Coil 


Figure  6.  Basic  schematic  of  type-II  NMR  quantum  computer.  The  gradient  coil  creates  a  gradient  in  the 
magnetic  field  so  that  the  nuclear  spin  energy  levels  are  shifted  by  different  amounts  depending  on  their 
physical  location  in  the  liquid  sample.  This  allows  the  RF  coil  to  address  different  parts  of  the  liquid 
sample  with  different  frequency  radio  pulses.  Each  group  of  molecules  that  the  RF  coil  can  address  with 
one  set  of  frequencies  is  a  node  in  a  type  II  quantum  computer.  In  each  node  there  are  many  molecules  that 
are  manipulated  simultaneously,  so  that  measuring  a  node  is  an  example  of  ensemble  measurement.  This 
figure  was  used  with  permission  from  [6]. 
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To  make  a  NMR  machine  a  type-II  quantum  computer,  the  liquid  sample  is 
effectively  split  into  nodes  using  a  magnetic  field  gradient.  This  gradient  splits  the 
nuclear  spin  states  by  different  amounts  depending  on  a  molecule’s  location  in  the  liquid 
sample,  allowing  one  to  address  different  sections  of  the  sample  with  different  frequency 
radio  pulses. 

The  FQLGA  has  two  properties  that  make  it  ideal  for  implementation  on  a  type  II 
NMR  quantum  computer.  First,  it  requires  no  more  than  three  entangled  qubits,  making 
it  possible  to  run  the  algorithm  on  a  type  II  computer  with  three  qubits  per  node.  NMR 
machines  with  three  qubits  per  node  have  already  been  successfully  demonstrated  using 
Alanine,  Trifluorobromoethylene,  and  Trichloroethylene  [15].  Secondly,  this  algorithm 
stores  information  in  the  probability  coefficients  of  the  basis  states.  This  information  can 
be  obtained  by  an  ensemble  measurement  over  all  the  molecules  in  a  node  of  a  NMR 
machine. 

3  Fluid  Dynamics 

Since  the  factorized  quantum  lattice  gas  algorithm  models  fluid  dynamics,  it  is 
worthwhile  to  briefly  review  this  subject  along  with  the  classical  lattice  gas  algorithms 
that  inspired  their  quantum  counterparts.  This  section  starts  with  a  brief  overview  of 
fluid  dynamics  before  discussing  classical  lattice  gas  and  lattice  Boltzmann  algorithms. 

3.1  Navier-Stokes  fluids:  macroscopic  scale 

The  following  section  follows  Landau  and  Lifshitz  [20],  Yepez  [4],  and  Buick 
[21],  The  long  wavelength  hydrodynamic  behavior  of  a  fluid  at  the  macroscopic  scale 
can  be  modeled  by  a  set  of  coupled  partial  differential  equations.  These  equations  model 
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mass  density  (p )  and  flow  velocity  ( u  )  fields,  and  are  called  the  continuity  and  Navier- 
Stokes  equations. 

Since  the  fluid  mass  change  in  a  region  51  comes  from  the  fluid  flux  through  the 
boundary  551 ,  p  and  u  must  obey 

dtp  +  di(pui)  =  0  (3.1) 

which  is  the  continuity  equation.  Here  the  shorthand  8t  =  d/8t  and  5.  =  5/5x;  is  used, 

along  with  Einstein  indicial  notation,  which  implies  summation  over  repeated  indices. 

The  field  equation  for  Newton’s  second  law,  which  expresses  the  change  in  the 
momentum  density  in  terms  of  the  stress  at  the  boundary  of  the  region  551 ,  is  Euler’s 
equation 

d,(put)  +  dj  n,=0  (3.2) 

where  the  momentum  flux  density  tensor  can  be  written 

n ij  =  P(p,  t)Sy  +  puiuj  -  cr '.. .  (3.3) 

The  first  two  terms  are  the  ideal  parts  of  the  momentum  flux  density  tensor,  which  are 
the  pressure  term  P(p,t )  and  the  convective  term  putt .  The  pressure  term  is  diagonal 
because  the  fluid  is  isotropic.  The  last  term  is  the  stress  tensor,  equal  to 
a  \j  =  q(8iuj  +  8jui  -  f  8kuk8tj )  +  < <f8ij8kuk ,  where  q  is  the  shear  viscosity,  (f  is  the  bulk 

viscosity,  and  D  is  the  number  of  spatial  dimension  of  the  system. 

Substituting  (3.3)  into  Euler’s  equation  gives  the  Navier-Stokes  equation 

p(dtui+uj8Jui)  =  -8iP  +  pvd2ui+  C  +  dfijUj  (3.4) 

v  cJ  J 
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where  v  =  yp  is  the  kinematic  viscosity.  This  equation  has  known  solutions  in  only  a  few 

simple  cases,  and  computer  modeling  with  various  numerical  techniques  are  typically 
necessary  to  solve  this  equation  for  more  complex  flows. 

The  kinematic  viscosity  v  is  a  measure  of  the  rate  of  decay  of  local  shears  in  a 
fluid,  and  determines  how  fast  a  fluid  will  relax  from  an  anisotropic  to  an  isotropic  flow 
field.  The  shear  viscosity  alone  is  responsible  for  the  damping  of  shear  waves  in  the 
momentum  density  field,  while  both  the  shear  and  bulk  viscosities  cause  damping  of 
compression  waves  in  the  mass  density  field. 

L  and  T  are  the  characteristic  length  and  time  scales  of  a  fluid  fluctuation. 
Examples  of  the  characteristic  length  for  a  hydrodynamic  flow  are  the  wavelength  of  a 
compression  wave  in  the  mass  density  field,  the  wavelength  of  a  shear  wave  in  the 
momentum  density  field,  or  the  diameter  of  a  vortex.  Examples  of  characteristic  times 
are  the  period  of  a  wave,  or  the  rotation  period  of  a  vortex.  The  mean  free  path  (A)  and 
time  (  t  )  are  the  average  distance  and  time  that  microscopic  particles  in  the  fluid  travel 
before  colliding.  Two  important  speeds  are  the  characteristic  flow  speed  v~L/T  and 
sound  speed  c  =  x/T . 

Relevant  dimensionless  numbers  are:  the  Knudsen  number  (Kn)  defined  as  the 
ratio  of  the  mean  free  path  to  the  characteristic  length,  the  Strouhal  number  (Sh)  defined 
as  the  ratio  of  the  mean  free  time  to  the  characteristic  time,  Mach  number  (M)  defined  as 
the  ratio  of  the  characteristic  velocity  and  sound  speed,  and  Reynolds  number  (Re) 
defined  as  the  product  of  the  characteristic  velocity  and  the  characteristic  length  divided 
by  the  kinematic  viscosity.  A  list  of  all  these  relevant  quantities  is  given  in  Table  4. 
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Table  4.  List  of  relevant  quantities  in  fluid  dynamics. 


Symbol 

Name 

Description 

P 

mass  density  field 

scalar  field  that  describes  the  fluid  mass  density 

u 

flow  velocity  field 

vector  field  that  describes  the  fluid  velocity 

t 1 

shear  viscosity 

causes  damping  of  compression  waves  in  mass  density 
field  and  shear  waves  in  momentum  density  field 

z 

bulk  viscosity 

causes  damping  of  compression  waves  in  mass  density 
field 

V 

kinematic  viscosity 

=  ri  /  p ,  determines  how  fast  perturbed  fluid  will  relax 

L 

characteristic  length 

length  of  fluid  perturbations 

T 

characteristic  time 

period  of  fluid  perturbations 

X 

mean  free  path 

particles’  average  distance  between  collisions 

T 

mean  free  time 

particles’  average  time  between  collisions 

V 

characteristic  speed 

~ LIT 

C 

sound  speed 

~  XI z 

Kn 

Rnudsen  number 

=  X/L 

Sh 

Strouhal  number 

=  t/T 

M 

Mach  number 

=  vl  c 

Re 

Reynolds  number 

=  oL!v  ~  M  /  Kn 

Returning  to  equation  (3.4),  one  can  see  that  the  one  dimensional  Navier-Stokes 
equation  simplifies  to 


du  du  d2u 

dt  dx  dx2 


(3.5) 


if  rj  =  £  =  P  =  0 .  This  is  a  simplified  model  of  turbulence  and  shock  formation  called  the 
Burgers  equation  [22],  In  section  4,  we  will  see  that  the  two  qubit  Factorized  Quantum 
Lattice  Gas  Algorithm  is  capable  of  accurately  modeling  this  equation. 


3.2  Classical  lattice  gas  algorithm:  microscopic  scale 

In  the  1980’s  a  class  of  algorithms  called  Lattice  Gas  Algorithms  (LG A)  were 
discovered  to  behave  like  a  Navier-Stokes  fluid  by  Wolfram  [23]  and  by  Frisch, 
Hasslacher,  and  Pomeau  [24],  raising  the  possibility  of  using  massively  parallel 
computers  running  LGAs  to  simulate  fluid  dynamics.  These  simulations  may  include 


21 


attractive  interactions  between  particles  to  create  multiphase  fluids  [25]  or  fixed  obstacles 
to  simulate  vortex  shedding  [26], 

Lattice  gas  algorithms  move  virtual  particles  along  an  imaginary  lattice  and 
change  the  particles’  momentums  using  collision  rules  when  they  meet  at  a  lattice  node. 
The  lattice  node  spacing  ( t )  is  defined  to  be  a  microscopic  scale  length  so  that  LGAs  are 
sometimes  said  to  model  fluids  at  this  scale.  In  fact,  lattice  gas  algorithms  grossly 
oversimplify  microscopic  particle  dynamics.  However,  this  turns  out  not  to  matter  since 
the  macroscopic  behavior  of  a  fluid  does  not  depend  directly  on  its  microscopic 
components.  This  is  evident  in  experiments  carried  out  using  wind  tunnels  and  water 
tanks  with  low  Mach  flows  and  similar  Reynolds  numbers,  since  the  results  of  both  types 
of  experiments  will  be  similar  [21].  Similarly,  LGAs  in  the  continuum  limit  (with  very 
small  i )  turn  out  to  be  accurate  models  of  fluid  dynamics. 

The  simulated  particles  in  a  lattice  gas  algorithm  are  located  on  the  nodes  of  a 
regular  lattice.  The  position  and  momentum  of  each  particle  is  specified  by  its  position 
on  the  lattice  and  a  displacement  vector.  The  displacement  vector  points  in  the  direction 
that  the  particle  will  move  at  the  beginning  of  a  time  step.  Particles  move  from  one  node 
to  another  in  a  process  called  streaming.  In  the  case  of  a  single  speed  lattice  gas,  all 
particles  move  at  the  same  velocity  c  =  \,  where  l  is  the  distance  between  lattice  sites 

and  r  is  the  time  step  interval  [23].  All  of  the  particles  stream  simultaneously  at  the 
beginning  of  a  time  step.  Most  LGAs  enforce  an  exclusion  principle  so  that  no  more  than 
one  particle  can  occupy  a  state  at  a  given  time,  though  there  is  usually  more  than  one 
particle  at  a  lattice  node.  A  state  is  defined  as  the  location  and  momentum  of  a  particle. 
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Each  state  is  typically  assigned  a  bit,  so  that  the  bit  of  a  full  particle  state  is  1  and  0  for  an 
empty  state. 

When  two  or  more  particles  meet  at  a  lattice  node,  their  momentums  change  in 
accordance  with  predetermined  collision  rules.  The  updated  particle  trajectories  are  then 
streamed  at  the  beginning  of  the  next  time  step.  This  process  is  shown  in  Figure  7,  with 
particle  trajectories  at  the  beginning  of  a  time  step  labeled  by  single  arrows  while  the 
trajectories  at  the  end  of  the  time  step  are  labeled  by  double  arrows. 


Figure  7.  Triangular  classical  lattice  gas  developed  by  Frisch,  Flasslascher,  and  Pomeau.  Particles  at  time  t 
are  marked  with  a  single  arrow;  those  at  the  next  time  step  t  +  T  are  marked  with  double  arrows.  Figure  is 
reproduced  from  [24], 

3.3  Classical  lattice  Boltzmann  algorithm:  mesoscopic  scale 

To  transition  from  the  microscopic  scale  Classical  Lattice  Gas  Models,  which 
contain  a  number  of  discrete  particles,  to  the  macroscopic  scale  Navier-Stokes  equation, 
which  contains  a  continuous  density  parameter  p,  it  is  necessary  to  convert  the  number  of 
particles  in  a  given  area  to  a  particle  density.  In  other  words,  particle  number  must  be 
replaced  by  a  continuous  statistical  particle  distribution  function.  This  is  analogous  to 
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describing  the  motion  of  a  group  of  microscopic  molecules  in  a  fluid  by  modeling  a 
mesoscopic  statistical  particle  distribution,  called  the  Boltzmann  distribution  f(x,u,t) , 
that  is  a  function  of  position,  velocity,  and  time. 

Boltzmann  mechanics  can  be  described  following  Gumett  and  Bhattacharjee  [27], 
where  we  consider  a  group  of  particles  /  (x,u,t)d:'xd3u  in  the  phase  space  volume 

element  d3xd3u  at  time  t.  These  particles’  positions  will  change  to  xr  =  x  +  ut  and  their 
velocities  to  u'  =u+^t  an  instant  later  at  time  t  +  z ,  so  that  they  occupy  a  new  volume 

in  phase  space:  d3x'd3u' .  Or  in  other  words, 

f(x  +  uT,u+^T,t  +  r)d3x'd3u'  =  f(x,u,t)d3xd3u.  (3.6) 

Any  change  in  the  particles  distribution  that  equation  (3.6)  does  not  account  for  must  be 
due  to  collisions  Q (/) ,  so  that1 

[/(x  +  ut,u  +  ^r,t  +  r)~  f(x,u,t)Jd3xd3u  =  Q.(f)d3xd3ur  .  (3.7) 

Expanding  this  result  in  a  Taylor  series  and  taking  the  limit  that  r  is  zero,  one  arrives  at 
the  well  known  Boltzmann  Equation 

^.+S.V/+-.V,/  =  n(/)  (3.8) 

ot  m 

In  contrast,  if  we  reconsider  the  discrete  space-time  of  a  lattice  gas  algorithm  and  set 
external  forces  equal  to  zero,  equation  (3.7)  becomes  the  finite  difference  equation 

fu(x  +  UT,t  +  T)-fu(x,t)  =  Q(fu)  (3.9) 


1  The  Jacobian  of  the  change  in  the  phase  space  volume  Jd3 xd3u  =  d3x'd3u'  is  equal  to  one. 
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where  we  now  choose  to  index  /  by  the  velocity  u  .  This  is  called  the  lattice  Boltzmann 
equation. 

There  are  two  methods  of  modeling  this  mesoscopic  equation.  The  first  approach 
is  to  directly  simulate  the  equation  on  a  lattice  using  continuous  values  for  the  particle 
occupation  of  a  state  instead  of  the  binary  1  for  “particle  present”  and  0  for  “no  particle.” 
Algorithms  that  follow  this  approach  are  called  lattice  Boltzmann  algorithms.  The 
second  approach  is  to  model  discrete  microscopic  particles  on  a  lattice  gas  simulation. 

The  governing  lattice  Boltzmann  equation  is  then  derived  as  an  approximate  description 
of  the  averaged  mesoscopic  dynamics  [21]. 

One  way  to  average  over  the  lattice  is  called  coarse  grain  averaging  and  works  by 
placing  a  mesoscopic  “superlattice”  over  the  microscopic  lattice  as  shown  in  Figure  8a, 
and  taking  the  average  occupation  probability  as  the  value  of/at  a  particular  superlattice 
site.  The  second  method  takes  an  ensemble  average  of  many  independent  microscopic 
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Figure  8.  a)  Coarse  grain  averaging  works  by  taking  the  average  over  all  the  microscopic  states  inside  the 
mesoscopic  superlattice,  b)  Ensemble  averaging  works  by  taking  the  average  over  many  independent 
microscopic  realizations  to  obtain  the  particle  distribution  at  each  site. 
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realizations  to  arrive  at  a  value  for  f  The  factorized  quantum  lattice  gas  algorithm  uses 
the  second  method  to  obtain  the  distribution / 

To  transition  from  a  discrete  (in  space)  mesoscopic  scale  simulation  to  a 
continuous  macroscopic  scale  Navier-Stokes  simulation  one  must  let  the  lattice  cell  size 
approach  zero  in  the  continuum  limit  as  shown  in  Figure  9.  In  this  limit,  it  is  possible  to 
perform  a  Chapman-Enskog  expansion  to  derive  the  macroscopic  governing  equation  of 
the  fluid  [24], 


Figure  9.  Decreasing  the  mesoscopic  lattice  cell  size  towards  zero  increases  the  simulation  resolution  and, 
in  the  continuum  limit,  will  approximate  a  continuous  macroscopic  field. 


When  performing  this  expansion,  one  should  note  that  the  particles  in  a  lattice  gas 
algorithm  undergo  random  walk.  That  is,  a  tagged  particle  will  move  a  distance  which 

asymptotically  approached L  =  ni^lhr  after  streaming  n2  times.  Since  particles  stream 
at  the  end  of  every  time  step,  n2r  =  T.  This  implies  that  in  random  walk  processes  Sh  ~ 
Kn2  ~  (1  In)2  =  s2  — a  condition  called  diffusive  ordering,  which  produces  viscous 
hydrodynamic  behavior  [5,  28],  This  will  become  important  as  the  Chapman-Enskog 
analysis  is  described  in  further  detail  in  the  next  section. 
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4  Factorized  Quantum  Lattice  Gas  Algorithm 


As  discussed  in  section  2.3,  type  II NMR  quantum  computers  are  well  suited  to 
simple  algorithms  that  requires  massive  parallelism.  LGAs  fit  this  description  well,  and 
it  was  this  observation  that  lead  Yepez  [2-5]  to  develop  the  factorized  quantum  lattice  gas 
algorithm  (FQLGA)  to  test  the  modeling  utility  of  quantum  computers.  The  algorithm  is 
called  “factorized”  because  it  is  not  meant  to  run  on  a  fully  coherent  computer,  but  rather 
on  one  that  is  made  up  of  many  smaller  quantum  computers.  Constant  measurement  of 
the  system  allows  one  to  transfer  classical  information  between  the  smaller  quantum 
computers  so  that  the  system  need  not  be  fully  coherent. 

Figure  10.  The  1-D  factorized  quantum  lattice  gas  model  developed  by  Yepez.  Each  lattice  site  is 
simulated  by  a  node  on  a  type  II  NMR  quantum  computer.  The  probability  of  finding  a  particle  moving 
right  at  lattice  site  l  is  given  by  qubit  1  in  node  l,  and  the  probability  of  finding  a  particle  moving  left  at  that 
lattice  site  is  given  by  qubit  2.  Since  there  are  many  computers  per  node  in  a  type  II  NMR  machine,  one 
can  perform  an  ensemble  measurement  on  each  node  to  obtain  the  probabilities  that  will  be  streamed  via 
classical  communications  channels  to  neighboring  nodes. 


The  FQLGA  is  similar  to  a  classical  LGA  in  that  it  simulates  particles  moving 
along  a  lattice  as  shown  in  Figure  10.  As  in  a  classical  LGA,  the  particles  move  via 
streaming  and  obey  collision  rules  when  two  particles  meet  at  a  node.  However,  unlike  a 
classical  LGA,  the  collision  rules  are  unitary  operations  which  mix  those  states  at  a  lattice 
site.  Following  the  collision  operation,  the  updated  probabilities  for  particles  moving 
right  (particle  1)  and  left  (particle  2)  are  obtained  via  ensemble  measurement.  These 
updated  probabilities  are  classically  streamed  to  nearest  neighbor  lattice  sites,  marking 
the  end  of  a  time  step.  Since  the  collision  operator  will  only  mix  those  states  at  a  given 
lattice  site,  each  lattice  site  can  be  simulated  by  a  node  on  a  type  II  quantum  computer. 
Thus,  the  1-D  FQLGA  developed  by  Yepez  with  two  particles  per  lattice  site  (moving 
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right  and  left)  requires  a  type  II  quantum  computer  with  two  qubits  per  node.  This  entire 
process  is  described  in  more  detail  in  the  following  section. 

4.1  The  four  steps  process  for  the  FQLGA 

4.1.1  Step  1:  computational  memory  state  encoding 

The  first  step  in  the  FQLGA  is  to  encode  the  computational  memory  state  for  the 
two  qubits  \ql(xl,t)')  and  \q2{xl,t)),  representing  particles  1  and  2  at  the  lattice  site  l 

(located  at  xl ),  at  time  t.  At  time  t0,  the  probabilities  are  given  by  initial  conditions 
provided  by  the  user.  Subsequent  probabilities  are  determined  by  the  algorithm.  The 
probability  for  particle  m  to  exist  at  a  lattice  site  located  at  xt  at  time  t  is  written 

pm  (x, ,  t) ,  so  that  the  mth  qubit  is  encoded  as 

kra(x/»o>=V^(x/’0|i>+Vi-^ra(^,o|o>.  (4.i) 

Here  the  basis  state  1 1)  means  a  particle  exists  in  the  simulated  state  and  1 0)  means  it 
doesn’t.  The  state  of  the  entire  node  is  called  the  local  state  and  is  given  by  the  tensor 
product  of  the  qubits  |  qx  (xl ,  t))  and  |  q2  (xl ,  t)) .  It  has  the  form 

\v(xi>  0)  =  1 4i  0/ ,  0)  ®  |  <h  (• xi  >  0) 

=  4PiPi  |11)  +  Va(1-/L)|10)  (4-2) 

+V(1“A)/?2  |oi)  +  V(i-a)(1-^)|oo) 

where  I  have  dropped  the  explicit  time  and  position  dependence  of  pm  (x, ,  t )  for 
convenience,  and  the  labels  inside  the  kets  are  ordered  for  particles  1  and  2.  Note  that  I 
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have  assumed  the  qubits  are  distinguishable  since  the  local  ket  is  neither  symmetric  nor 
antisymmetric1. 

4.1.2  Step  2:  collision 

Following  memory  state  encoding  the  local  state  undergoes  unitary  evolution. 

This  is  analogous  to  the  collision  operator  in  the  classical  LGA  and  is  therefore  labeled 

C .  Thus,  the  local  state  becomes 

\vXx„t))  =  C\V(x„t)).  (4.3) 

If  we  use  the  basis  in  (2.2)  and  choose  a  collision  operator  that  conserves  particle  number 
then  the  operator  will  be  block  diagonal  and  have  the  form 

^10  0 

0  e'V*  cos  0  e'Vy'  sin  6 
0  -elGe~ul>  sin 6  ea e  cos 6 

v0  0  0 

The  block  that  mixes  the  states  |l0)  and  |0l)  is  a  U(2)  matrix. 

4.1.3  Step  3:  measurement 

This  step  destroys  the  quantum  superposition  and  measures  the  probability  of 
each  qubit  to  be  in  the  state  |l) .  As  discussed  in  section  2.1.2,  this  probability  can  be 
obtained  via  an  ensemble  measurement,  and  for  convenience  may  be  expressed  as 

P  'm  ( xi >■ 0  =  {v  Xxi » ■ 0 1 K  |  yr  \x, ,  0)  (4.5) 

1  This  is  reasonable  since  one  designs  a  quantum  computer  so  that  its  qubits  are  distinguishable.  For 
instance,  the  nucleons  chosen  to  represent  qubits  have  different  spin  energy  levels  in  a  NMR  computer. 


0^ 

0 

0 

1 


(4.4) 


29 


where  the  operators  nm  are  given  in  (2.8)  and  (2.9).  Notice  that  the  equation 

Pm  ( xi  >  0  =  {wix,  ,t)\nm\i/f(xlf  t )>  (4.6) 

also  holds.  This  is  of  no  consequence  now  but  will  be  useful  later. 

4.1.4  Step  4:  streaming 

Each  lattice  site  is  updated  following  steps  1  through  3  and  the  resulting 
probabilities  are  streamed  to  the  adjacent  lattice  sites  via  classical  communications 
channels  so  that 

Pm ( xi  +ej,  t  +  T)  =  p'm (x, ,  t )  (4.7) 

where  e\  =  1  for  the  right  streaming  qubit,  ei  =  -\  for  the  left  streaming  qubit,  and  t  is  the 
lattice  spacing.  This  signals  the  end  of  a  time  step  r,  after  which  the  entire  process  is 
repeated.  The  simulation  will  include  many  time  steps  and  is  typically  completed  when 
the  system  reaches  equilibrium  and  exhibits  no  further  change. 

4.2  Quantum  lattice  Boltzmann  equation 

By  simple  substitution,  all  four  steps  in  the  FQLGA  can  be  encapsulated  in  one 
equation.  This  is  carried  out  as  follows: 

Pm(.xi+emt>  t  +  T)  =  P'm(xn  0 

=  {vX*i>t)\hm\y/'{xl,f))  (4.8) 

=  {ip(X/,  0|CfnmC|^(A>  0). 

With  a  few  modifications,  this  equation  becomes  the  finite  difference  quantum 
lattice  Boltzmann  equation,  analogous  to  the  classical  equation  (3.9).  The  first  step  is  to 
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reinterpret  the  probabilities  pm (x, ,  t)  as  a  mesoscopic  Boltzmann  field  pm (x, ,  t )  = 

■  Therefore,  equation  (4.6)  can  be  rewritten 

fm ( xi > 0  =  {y(x, ,  0 1 nm  \t//(x, , 0)  (4.9) 

and  the  result  of  (4.8)  can  be  rewritten 

fm(Xl+emt,  t  +  T)  =  (i//(xl,  t)\CfnmC\i//(x„  t))  .  (4.10) 

Then  all  that  is  left  to  do  is  to  subtract  equation  (4.9)  from  (4.10)  so  we  obtain 

fm(xi+eJ >  t  +  T)-fm(xi’  t)  =  (Y(x„  0|(C’t«raC-nm)|^(x;,  0)  n) 

This  is  the  quantum  lattice  Boltzmann  equation.  This  can  be  further  expanded  by 
inserting  the  vectors  for  {ip(xn  t)\  and  |^/(x/,  t))  along  with  the  matrices  for  C  and  nm 
into  (4.1 1).  With  much  algebraic  manipulation  [12],  the  collision  function  becomes 

nu=T(sinJ«[(l-/2)/;-(l-y;)/2]  +  sin(26i)cos(^-^(l-/i)/2(l-/2))(4.12) 
or,  written  more  simply 

fi1>2  =  +  (sin!  ew,  -  A )  +  sin(2$)  cos(,d  -  ^/(l  ■ -  /  )/2  (1  ■ -  /2 ) ) .  (4.13) 

4.3  Chapman-Enskog  expansion 

As  mentioned  in  section  3.3,  one  can  derive  the  macroscopic  governing  equation 
in  the  continuum  limit  of  the  lattice  by  performing  a  Chapman-Enskog  expansion  of  the 
lattice  Boltzmann  equation.  This  section  will  follow  Yepez  [4,  10]  to  derive  the 
macroscopic  equation  for  his  model. 
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The  Chapman- Enskog  expansion  works  by  taking  a  Taylor  series  expansion  of  the 
lattice  Boltzmann  equation  around  local  equilibrium.  In  physical  systems,  local 
equilibrium  is  the  state  where  particles  are  in  thermodynamic  equilibrium  with  one 
another  across  mesoscopic  or  microscopic  scale  lengths.  In  classical  lattice  Boltzmann 
and  quantum  lattice  gas  algorithms,  local  equilibrium  is  obtained  at  a  lattice  site  when  the 
collision  function  no  longer  changes  the  particle  distribution  at  that  site. 

We  expand  around  local  equilibrium  because  at  the  mesoscopic  scale,  most 
systems  are  at,  or  very  near  to,  thermodynamic  equilibrium.  It  is  only  at  the  macroscopic 
scale  that  there  are  free  thermodynamic  variables  such  as  local  density,  temperature,  and 
momentum.  Thus,  a  macroscopic  description  of  a  fluid  comes  from  a  patchwork  of 
slowly  varying  systems  at  or  very  near  equilibrium  [24], 

The  dimensionless  numbers  Kn  (Knudsen  number),  Sh  (Strouhal  number),  and  Re 
(Reynolds  number)  discussed  in  section  3.1  can  be  used  to  determine  how  close  a  system 
is  to  equilibrium — at  equilibrium  these  numbers  vanish.  For  instance,  at  equilibrium  the 
characteristic  length  scale  is  infinitely  large  compared  to  mean  free  path,  so  Kn  ~  0. 
However,  hydrodynamic  behavior  is  also  attained  in  the  long  wavelength  limit  where 
these  numbers  are  close  to  zero.  Thus,  it  should  be  of  no  surprise  that  expanding  the 
lattice  Boltzmann  equation  around  local  equilibrium  should  result  in  hydrodynamic 
behavior. 

In  the  following  section  the  local  equilibrium  value  of fm  is  derived.  In  section 
4.3.2  this  result  is  used  in  the  Chapman-Enskog  expansion  of  the  quantum  lattice 
Boltzmann  equation. 
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4.3.1  Local  equilibrium 


For  simplicity,  we  will  label  the  equilibrium  values  of fm  as  dm.  Local  equilibrium 
is  defined  as  the  condition  where  dm (x;  +em£,  t  +  r)-dm(xl,  t)  =  0;  that  is,  when 
collisions  cause  no  further  change  in  the  distribution  function^.  From  this  we  can  see 
that  Q„,  | |/|  2=f/|  2  =  0 ,  or  from  (4. 12) 

sin2  6\{\-  d2)d{  - (1  -  dx )d2 ]  +  sin(2 0) cos(^ -  (1  - ~dx )d2 (1  —  d2)  =  0  (4.14) 

Dividing  this  equation  by  (1  -  Jj)(l  -  d2 )  and  rearranging  we  get 


dx 


(1  ~d2)  (1  ~dx) 


=2  cot(0)  cos(^  -  £)  '  ^ 


n-d^d-d,) 


(4.15) 


Taking  the  equilibrium  probabilities  to  have  the  form 


dx  = 


yz  + 1 


and 


d2  =  ■ 


+  1 


(4.16) 


and  substituting  this  into  (4.15)  with  some  manipulation  gives  the  quadratic  equations 

y2  +  lay  -1  =  0  (4.17) 

where  a  =  cot  0  cos(^  -%).  We  take  the  positive  root  solution  of  this  (so  that  dm  will  be 


positive)  so  that 


y  =  yfa2  +1  +a 
—  =  \l~cc2  +1  -  a. 

r 


(4.18) 


Noting  that  the  total  number  density  at  a  lattice  site  is  dx+d2=p  and  substituting  (4. 16) 
and  (4.17)  into  this  expression,  we  obtain  a  quadratic  equation  in  z 


pz  + 


1 

y  +  — 

V  rj 


(p-l)z  +  (p-2)  =  0. 


(4.19) 
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When  the  positive  root  solution  of  (4.19)  along  with  (4.17)  is  substituted  into  (4.16)  one 
finally  arrives  at 


du2  =f  +  ^(VW-Vl  +  «2(P-1)2)  (4.20) 

which  are  the  equilibrium  values  for  the  two  qubit  quantum  lattice  Boltzmann  equation. 

4.3.2  Taylor  series  expansion  around  local  equilibrium 

To  keep  track  of  the  order  of  the  expansion,  one  uses  a  “smallness  parameter”  s . 
This  is  defined  to  be  on  the  order  of  the  Knudsen  number:  Kn  =  i  /  L  ~  e .  Like  the 
classical  lattice  gas  algorithm,  the  factorized  quantum  lattice  gas  algorithm  obeys 
diffusive  ordering.  Therefore,  s2  must  be  on  the  order  of  the  Strouhal  number  Sh  = 
f~s2  [4],  The  QLBE  can  thus  be  written 

fm(X[+£eJ,  t  +  £2t)~  fm(xn  t)  =  Clm.  (4.21) 

With  this  we  can  now  find  the  Taylor  series  of  the  left  side  of  (4.21),  which  is 
just  the  quantum  lattice  Boltzmann  equation  rewritten  to  include  £ .  This  gives 

£2T^L  +  £ej^L+ e2e2J2  - +  0(£3 )  =  O  .  (4.22) 

dt  m  dx  m  2  dx2  K  ’  m  y  ’ 

We  did  not  expand  the  right  hand  side  of  (4.21)  because  it  turns  out  to  be  easier  in  the  end 
not  to.  Notice  that  the  equation  is  still  exact,  since  we  implicitly  keep  higher  order  terms 
in  the  factor  0(£-3)  . 

The  crucial  ansatz  is  to  now  assume  that  fm  can  be  expanded  around  local 
equilibrium  as  the  sum  of  dm  and  a  small  deviation  3  fm ,  which  may  in  turn  be  expanded 
in  powers  of  £ : 
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f.-dm+£Sfr+s2Sf«'+  0(£J).  (4.23) 

Inserting  this  into  (4.22)  and  explicitly  writing  out  only  those  terms  up  to  second  order  in 
s  we  obtain 


ijdm .  __  „a*„  .  „2_  nesC  .  _2„2„2 1  s-dm . 


£  T 


~^!L  +  £eJ^L  +  £zeJ — ^ 
St  m  dx  m  dx 


-  +  £  <?,„ 


2  dx2 


+  0(0  =  0, 


(4.24) 


This  is  the  QLBE  expanded  around  equilibrium.  However,  we  do  not  have  an 
expression  for  8 /^0) .  To  obtain  one,  we  take  the  first  moment  of  the  QLBE  and  solve  for 


8 /j0) .  First,  expand 


Q  =  Q 

m  m 


dQ„ 


^  dfn 


8C+  0(£z) 


(4.25) 


f\,2~d\. 


Note  that  the  first  term  in  this  expression  is  equal  to  zero  from  the  definition  of  local 
equilibrium.  Equating  like  powers  of  £ ,  the  first  moment  of  (4.21)  is  therefore 


.  8dm  dQ„ 

£eJ^  =  £YJ~ 


dx 


r  dfn 


Sfn 


(0) 


(4.26) 


f\ ,2  ~d\ . 


From  equation  (4.13)  we  see  that  Qj  =  -Q2  =  Q .  Then,  for  simplicity,  equation 
(4.26)  can  be  rewritten  in  vector  form 


ele 


dcT 
dx  i 


=  £j\8f 


(0) 


(4.27) 


where  the  vectors  ^ 


and  1 8 /(0) 


are 


dd,r 


dx 


(to 8\ 

dx 

dd2 

V  dx  J 


and 


\8f(0)}  = 


(4.28) 


and  the  matrices  J  (the  Jacobian)  and  e  are 
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(  6Q 

df\ 

8Q  \ 

df2 

r  ^ 

V 

-dQ 

V  dfi 

-dQ 

¥2  ) 

II 

•23" 

11 

~^2  y 

and 


f  1 
,0 


0^ 

-h 


(4.29) 


Solving  for  1 8 /(0)  j  is  not  as  easy  as  finding  the  inverse  of  J ,  however,  because 
this  matrix  is  singular.  Yepez  utilizes  two  (equivalent)  methods  to  find  a  consistent 
\Sf(0))  [10].  The  first  is  to  multiply  both  sides  of  (4.27)  by  J .  The  second  method  is  to 

multiply  both  sides  of  (4.27)  by  a  “generalized  inverse”  J~ ,  which  Yepez  has  invented. 
This  matrix  is  similar  to  the  Moore-Penrose  pseudoinverse  [29], 

The  eigenvalues  and  left  and  right  eigenvectors  of  J  are 


4=0  (£,|  =  (1  1) 


1  \  1 

"y " 

Ei)  = 

J  -J 
u\  u  2 

\J\) 

I  J 2 


E2  = 


J2  Jx 


(y  j2) 


e2  = 


f  i  A 


v-ly 


(4.30) 


where  the  right  eigenvectors  (often  simply  called  eigenvectors)  satisfy  J  |  Em  'j  =  Xm  |  Em  ^ , 
the  left  eigenvectors  satisfy  (Em  |  J  =  Xm  (Em  | ,  and  the  eigenvector  lengths  are  selected  so 


that  (E\E\  =  8,.,„ .  From  this  we  see  that  J  is 


m  I  n  mn 


J  =  A,\Ei)(Et\  +  A1\E2)(E2\  =  A1\E1)(E2\  = 


f  J  J  ^ 


-y  -y 


(4.31) 


2  J 


since  \  =  0 .  This  equation  is  analogous  the  spectral  decomposition  of  a  Hermitian 
matrix,  and  is  equivalent  to  square  matrix  diagonalization.  Any  n  by  n  square  matrix  M 
with  n  independent  eigenvectors  can  be  diagonalized  using  S 'MS  =  A  [30],  The 
columns  of  S  are  the  right  eigenvectors  of  M,  A  is  a  diagonal  matrix  with  corresponding 
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eigenvalues,  and  the  rows  of  S  1  are  the  left  eigenvectors  of  M  since  (Em  |  En )  =  8mn . 

Therefore  M  =  SAS  1  is  equivalent  to  (4.31). 

From  (4.31),  it  is  clear  that 

J2=%\E2)(E2\  =  A2J.  (4.32) 

As  was  mentioned,  the  first  method  Yepez  uses  to  find  8 is  to  multiply  both  sides 


of  (4.27)  by  J .  Then  one  obtains 

Jle  yj  =  J2\8f{0)) 

Jle  ^  =  JA2\8f(0) 

which  implies 


(4.33) 


The  second  (and  equivalent)  method  Yepez  uses  is  to  multiply  both  sides  of  (4.27) 
by  his  generalized  inverse.  The  generalized  inverse  is  analogous  to  the  inverse  of  a 
nonsingular  square  matrix  M”1  =  SA_1S_1 .  Yepez  uses  an  identical  construction  for  his 
generalized  inverse  except  that  he  replaces  A"1  with  a  matrix,  A“*n ,  in  which  only  the 


nonzero  diagonal  components  are  inverted.  The  procedure  for  constructing  the  Moore- 
Penrose  pseudoinverse  is  similar,  except  the  vectors  which  make  up  S”1  and  S  are  the  left 


and  right  eigenvectors  of  MMf .  It  can  be  shown  that  when  M  is  invertible,  the  least 
squares  solution  for  Mx  =  b  is  *  =  M;L/0b,  where  M  psuedo  's  the  Moore -Penrose 


pseudoinverse  [30]. 
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Yepez’s  generalized  inverse  for  J  is 


a2  a,  a 


(4.34) 


Multiplying  both  sides  of  (4.27)  from  the  left  by  this  generalized  inverse,  one  obtains 


/  1  te 

J  gen1 


\jie 


Jle 


dd 

dx 

dd 

dx 

dd 

dx 


~  =  J2\5fw 


(0) 


(4.35) 


which  may  be  further  simplified  following  the  same  steps  given  in  (4.33).  Thus,  the  two 
methods  are  equivalent.  If  one  takes  the  solution  obtained  in  (4.33)  and  substitutes  it  into 
(4.24)  one  arrives  at 

d2d„ 


,2_ddm  ,  _  „ddm  ,  _2  1  _2„2  d %  ddn 


s  r — -  +  seml- 
dt 


+  £z^ezJz 
dx  An  dx  dx 


+ 


f\  1N 

- 1 - 

A2  2 


dx2 


-  +  ()(£•)  =  Qm  (4.36) 


The  next  step  is  to  sum  these  equations  over  m,  noting  that  d]  +d2  =  p ,  and 


e2m  =1.  This  gives 


,2 ,dp  ,  ,  d  2  £  dl1  dp 


s  t  —  +  s£  —  (dx  -d2)  +  s 


2^7  3  ~  f  J  j  ^ 


dt 


dx 


dx  dx 


+ 


\^2 


£2i2^-  +  O(£3)  =  0 


dxz 


(4.37) 


In  what  follows,  Yepez  restricts  himself  to  small  a  to  simplify  the  resulting  equations. 
That  is,  he  assumes  that  the  angles  6 ,  <j) ,  and  £  in  the  collision  operator  are  such  that 

\a\  =  cot#cos(^  -  £)|  «  1 .  Then,  using  the  equilibrium  values  (4.20),  the  second  term  in 


(4.37)  is 


38 


(4.38) 


s£—(dl-d2)  =  s£—  — -(Vl  +  a2  -Ji  +  a2(p-l)2 
dx  dx  a  V 


=  -s£a(p-\)^-(\  +  a\p-\)2)  1/2 
=  ^a(l-p)|^(l  +  0(a2)) 


/  2  2 \— ^  ^ 

where  the  Taylor  series  expansion  of  ( 1  +  a  (p  - 1)  1  with  respect  to  a  was  taken  in 


the  last  line.  Calculating  the  components  of/ one  obtains 


f 


Jl2  =  sin2  0 


(2/j  j  l)/2  l(l  ^2,1 ) 
+1  -a—j=+  = = 

v  ^^(l-Zj)/^!-^) 


(4.39) 


This  implies  that 

A2=Jx-J2=-2  sin2  #(1  +  a2/ ( a ,  p))  (4.40) 

where  / (a,/?)  is  very  complicated  but  has  the  important  property  that 
f(a,p)  =  1  +  0(a) .  The  resulting  equation  is 

^ +  -  cot  0  cos (<f>  -  £)(1  -  p)  ^  =  cot2  6» — ^4  +  0(f 3 ,  «z2 )  (4.41) 

dt  r  dx  t  dx 

Dropping  the  terms  implicit  in  0(s3,sa2)  one  obtains  a  partial  differential  equation  that 

models  the  FQLGA  in  the  continuum  limit  of  the  quantum  lattice  Boltzmann  equation 

(4.1 1),  accurate  to  first  order  in  time  and  second  order  in  space  for  small  a  . 

For  u=cs(\- p)  where  cs  =  |-cot#cos(^-£)  is  the  speed  of  sound  and 


v  =  ycot2  0-y  is  the  kinematic  viscosity,  equation  (4.41)  becomes  the  Burger’s  Equation 
introduced  in  section  3.1. 


du  du  d2u 
—  +  u  —  =  v — - 
dt  dx  dx2 


(4.42) 
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Thus  we  have  completed  the  Chapman-Enskog  expansion  by  taking  the  Taylor  series 
expansion  of  the  quantum  lattice  Boltzmann  equation  around  local  equilibrium,  valid 
when  |cot  <9cos(^  -  £)| « 1 .  This  entire  process  can  be  a  bit  difficult  to  follow  so  it  is 
summarized  here: 

1 .  Expand  fm  around  local  equilibrium:  fm=dm+  sS /J;0)  +  ... 

2.  Insert  the  expanded  fm  into  a  Taylor  series  expansion  of  the  QLBE,  explicitly  writing 

out  only  those  terms  of  order  s 2  or  lower  (since  the  first  derivative  in  time  is  on  the 
order  of  s 2 ). 

3.  Use  the  first  moment  of  the  QLBE  to  solve  for  the  unknown  S in  terms  of  dm  . 

Placing  the  equations  in  matrix  form  can  help,  but  it  is  nevertheless  tricky  since  the 
matrix  J  will  be  singular. 

4.  Insert  8 /^0)  into  the  results  of  step  2. 

5.  Sum  the  results  of  step  5  over  m  and  simplify,  taking  advantage  of  the  fact  that 

p=IX  • 

m 

The  resulting  equation  will  be  the  governing  partial  differential  of  the  quantum  lattice 
Boltzmann  equation  in  the  continuum  limit. 

4.4  Numerical  and  experimental  simulation  of  the  Burgers  equation 

Yepez  has  run  a  numeric  simulation  of  this  algorithm  with  a  lattice  length  equal  to 
256f  on  a  conventional  desktop  computer,  and  compared  its  results  to  a  known  analytic 
solution  of  the  Burgers  equation  [10].  The  initial  condition  of  the  simulation  was  a 
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sinusoidal  wave,  which  generates  a  shock  front  at  later  times.  For  the  simulation  he 
chose  r  =  1 ,  £  =  1,  0  =  n  /  4 ,  and  <j>  =  £  . 


t  32 
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Figure  11.  Numerical  results  of  the  Burgers  equation  simulation  carried  out  by  Yepez,  along  with  the 
analytic  solution.  Agreement  between  the  simulation  and  the  analytic  solution  is  generally  very  good,  with 
slight  deviations  occurring  at  later  times  as  a  steep  shock  front  is  formed.  The  simulation  shock  front  is 
sharper  than  the  analytic  shock.  This  figure  was  used  with  permission  from  [10]. 
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Figure  12.  Experimental  results  of  the  two  qubit  FQLGA  simulating  the  Burgers  equation  carried  out  on  a 
type  II NMR  quantum  computer.  The  black  dots  are  the  experimental  results  and  the  solid  gray  line  is  the 
analytic  solution.  This  figure  was  used  with  permission  from  [1 1], 


The  results  of  this  simulation  are  presented  in  Figure  11.  The  agreement  between 
the  simulation  and  analytic  solution  is  generally  very  good,  although  there  is  some 
divergence  at  later  times  as  the  steep  shock  front  forms.  The  simulation  produces  a 
sharper  edged  shock  than  the  curved  edge  analytic  solution.  The  agreement  is 
nevertheless  impressive  since  the  shock  front  appears  to  be  greatly  under-resolved. 
Comparable  classical  algorithms  used  to  model  the  Burgers  equation  require  significantly 
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more  lattice  sites,  216  f  =  655361 ,  and  time  steps,  218r  =  262 144 r ,  to  model  a  shock 
formation  with  this  accuracy. 

This  simulation  has  also  been  run  on  a  working  two  qubit  per  node  type  II NMR 
quantum  computer  with  16  nodes.  The  results  of  this  are  shown  in  Figure  12.  The 
dominant  errors  in  this  simulation  come  from  errors  in  applying  the  collision  operator, 
which  accumulate  over  time. 

5  Three  Qubit  FQLGA  Using  Most  General  Collision  Matrix 

The  two  qubit  FQLGA  developed  by  Yepez  can  be  extended  in  one  dimension  by 
adding  on  an  additional  qubit  per  lattice  site  that  does  not  move  during  streaming.  Thus, 
this  sort  of  algorithm  will  have  three  particles  per  lattice  site:  particle  one  streams  right, 
particle  two  does  not  move,  and  particle  three  streams  left.  The  difficulty  in  extending 
this  algorithm,  however,  is  that  each  additional  qubit  greatly  increases  the  complexity  of 
the  most  general  collision  operator  that  conserves  particle  number.  This  makes  the 
analytic  treatment  increasingly  difficult  to  carry  out.  Nevertheless,  I  have  derived  the 
quantum  lattice  Boltzmann  equation  for  the  most  general  three  qubit  collision  operator 
that  conserves  particle  number,  as  well  as  the  diffusion  equation  for  a  specific  collision 
matrix. 

In  what  follows,  I  will  discuss  the  analytic  treatment  of  the  most  general  three 
qubit  operator  before  deriving  the  diffusion  equation  in  full  detail.  Then  I  will  show  the 
results  of  numeric  simulations  of  the  three  qubit  FQLGA  using  the  collision  matrix  which 
models  the  diffusion  equation.  I  will  compare  this  simulation  to  the  analytic  solution  and 
investigate  the  convergence  and  numerical  stability  of  the  simulation. 


43 


5.1  Microscopic  scale:  matrices  and  basis  states 

Like  the  two  qubit  algorithm  developed  by  Yepez,  my  three  qubit  algorithm  is 
designed  to  conserve  particle  number.  For  this  reason,  there  may  only  be  mixing  between 
the  states  1 1 1 0) ,  |lOl) ,  and  1 0 1 1) ,  and  separately  1 00 1) ,  1 0 1 0) ,  and  jlOO) .  I  therefore 
choose  the  basis  given  in  (2.5),  so  that 

(5.1) 

where  pm  is  of  course  the  probability  of  finding  a  particle  in  state  m.  This  choice  makes 

the  collision  matrix  block  diagonal.  Thus  it  must  have  the  form 

10  0  0 
0  U(  3)  0  0 

0  0  U(3)  0 

0  0  0  1 

This  is  shorthand,  since  C  must  be  an  eight  by  eight  matrix.  The  entries  U(3)  represent 
the  most  general  three  by  three  unitary  matrices.  In  fact,  this  collision  matrix  can  be 
simplified  even  further  since  U(3)  =e‘a  SU(3).  From  equation  (4.10),  we  see  that  the 

/V  ,  /V 

system  dynamics  are  determined  by  the  matrix  C'nmC ,  where  nm  gives  the  probability  of 
finding  a  particle  in  state  m  and  represents 
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/v  .  /V 

Since  hm  is  diagonal  we  know  that  C'nmC  must  have  the  form 


f  1  0 

t  „~ia  . 


CfnC  = 


0  SU(3ye-,lT  -Dmy  emSU(3) 


0 

0 


0 

0 


SU(3jte-kr-Dm2-e!<rSU(3)  0 


0^ 

0 


0 


0 


(5.4) 


where  Dm  n  are  diagonal  matrices  that  depend  on  nm .  For  instance,  for  m  =  1 


fl 

0 

0^ 

^0 

0 

0^ 

0 

1 

0 

and  Dl2  = 

0 

0 

0 

0 

0 

OJ 

,0 

0 

1, 

(5.5) 


Obviously,  the  terms  e  m  can  be  factored  through  the  D  matrices  and  cancel  the  terms 


eia ,  so  without  any  loss  of  generality  one  can  use  the  simplified  collision  matrix 
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(5.6) 


'1  0  0  0^ 

0  SU( 3)  0  0 

0  0  SU(  3)  0 

v0  0  0  1, 

The  matrix  SU(3)  has  eight  free  parameters  and  is  too  complicated  to  write  down  here  but 
is  presented  in  Appendix  A. 


5.2  Mesoscopic  scale:  quantum  lattice  Boltzmann  equation 

The  quantum  lattice  Boltzmann  equation 

fm(xi+em t  +  T)-fm(x„  t)  =  (y/(x„  t)|(ctnmC-nm)|^(x;,  t)}  (5.7) 


developed  in  4.2  is  still  valid  for  my  three  qubit  algorithm.  Due  to  the  complexity  of  the 
most  general  SU(3)  matrix,  it  is  useful  to  temporarily  replace  it  with  the  over¬ 
parameterized  matrix 


r aeld“ 

be0" 

ce>9‘  ^ 

SU(3)  -> 

dew“ 

feWf 

iOg 

ge 

(5.8) 

he'0" 

V 

j*1 

keWt  J 

where  I  have  written  the  entries  of  the  SU(3)  matrix  as  complex  numbers  in  polar  form, 
parameterized  by  the  real  numbers  a-d,f-h,j,k  and  0u_d  f_h  jk .  Since  the  rows  and 


columns  of  a  unitary  matrix  are  orthonormal,  the  constraints 


1)  a2  +b2  +c2  =  1 

2)  d2+f2+g2=  1 

3)  h2+j2+k2=  1 

4)  a2+d2  +  h2=  1 

5)  b2+f2+j2=  1 

6)  c2+g2+k2=  1 
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7)  adei(e°~dd)  +  bfe(db~df)  +  cge =  0 

8)  ahei('9a~6^  +  bjei(db~dj)  +  ckei(9c~9k)  =  0 

9)  hdei(9"-9d)  +  jfe 

1 0)  abei{6“~6h)  +  dfe'Ud~Uf'  +  hje 

11)  acel('9‘,-9c)  +  dge 

12)  bcei{9b~9c)  +  fgeK°f~°s)  +  jke  drdt)  =  0 


+  kgem~9*)  =  0 
+  hjei(9h~9i)  =  0 
■'(*.-*«)  JU  Jnjw*-9*')  +  hkei(9h~9k)  =  0 


(5.9) 


must  be  true.  Note  that  the  conjugates  of  identities  7  through  12  must  also  be  true. 
Though  rewriting  the  SU(3)  matrix  in  the  form  (5.8)  appears  to  have  complicated  matters, 

/V  .  /V 

it  significantly  simplifies  the  matrix  multiplication  ( C'nmC - hm ) .  For  m=  l,  this  matrix 
is  equal  to 
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.(5.10) 


Then,  using  the  identities  4  through  6  and  10  through  12  in  (5.9),  this  matrix  becomes 


>■«%-*,) 
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0  -hkem~°k) 
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-‘(.o.-e,) 
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hkei{0k~0k)  jkei(0J~0k) 
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jke~K0J~0k)  0 


-1  +  k2 
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(5.11) 


Multiplying  this  matrix  from  the  left  by  (y  |  and  from  the  right  by  |  if/')  and  using  identity 
3  from  (5.9)  one  obtains  the  collision  function  for  the  first  QLBE  collision  function. 


47 


(5.12) 


Qi  -  ( k 2  -1  )pt  +  j2p2  +  h2p3  + 

2hj cos \0h  - 0j](  1  - 2 px )yj(l  -p2)p2(l-p3)p3  + 

2 hk cos[0h  - 0k ](1  - 2 p2 )^(1- a)a(1~  A)A  + 

2 jk cos [6^.  -0k](\- 2p3 )yl(l-pl)pl(l-p2)p2  ■ 

With  similar  work  one  obtains  the  second  and  third  collision  functions: 


Q2  -  g2  px  +  (/2  - 1  )p2  +  d2p3  + 

2df  cos[0d  2pl )yjt (1  - /?2 )/?2 (1  - ~p2 )/?3  + 

2 dg cos[^  - 6g ](1  - 2jp2)V(1-A)A(l-ft)/?3  + 
2/gcos^  -  <9g  ](1  -  2^  y  (1  - ~p{  XA  (1  -  A )  A 

=  c  p^-\-  b  P2  +  ($  —  1) p^  + 

2 ab  cos[0a  -  9b  ](1  -  2/?,  )V(l-/?2)/?2(l-/73)/?3  + 

2ac cos[(9fl  -  0C  ](1  -  2p2)pj  -  a  ) A  (1 -  A ) A  + 

26c cos[6»,  -  0C ](1  - 2/?3 )V(1-a)a(1-A)A  • 

Using  identities  10  through  12  in  (5.9)  and  their  conjugates  one  can  see  that 

13)  ab  cos (0a  -0b)  +  df  cos(0d  -  0f )  +  hj  cos (0h  -0j)  =  O 

1 4)  ac cos(0a  -0c)  +  dg cos(0d  -0g)  +  hk cos {0h  ~0k)  =  0 

15)  be  cos  (0b  -0C)  +  fg  cos(<9/  -  0g )  +  jk  cos {0.  -0k)  =  O 


(5.13) 


(5.14) 


so  that  Qj  +  Q2  +  Q3  =  0 . 


5.3  Macroscopic  scale:  Chapman-Enskog 


As  discussed  in  section  4.3.2,  the  quantum  lattice  Boltzmann  equation  can  be 
expanded  in  a  Taylor  series  around  local  equilibrium  to  yield 


s2z 


dd 

_ m 

dt 


-  +  se„ 


.  dd,n 

dx 


■  +  £  el 


ssp 

dx 


(0) 


-  +  £  e. 


2„2  „2 


2  dx2 


+  0(*3)  =  Qm. 


(5.15) 
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One  may  sum  these  equations  over  m  with  ex  =  1 ,  e2  =  0,  e3  =  - 1 ,  and  p  =  dx+  d2+  d3  to 
obtain 


e2T^-  +  st-f(dl-d,)  +  e2t-^(Sflm-Sf,m)  +  s2e2\^-  +  O(e,)  =  0.  (5.16) 
ot  ox  ox  2  ox 

The  key  difficulty  now  is  finding  the  local  equilibrium  values  dm ,  which  one 


needs  to  solve  for  (dx  -  d3 )  and  (Sfx0)  -8 /3<0))  to  complete  the  Chapman-Enskog 


expansion.  Since  the  collision  functions  sum  to  zero,  only  two  of  these  functions  are 
linearly  independent.  The  collision  functions  are  each  equal  to  zero  at  local  equilibrium. 
Therefore,  the  three  equations  that  one  must  solve  to  obtain  the  three  unknown 
equilibrium  values  are 


Q,l  , 

=  0 

Qj 

=  0 

2  |  Pm  =dm 

p  —  t/j  +  d2  +  d2  . 


(5.17) 


The  complexity  of  the  collision  functions  necessitates  making  a  variable 
substitution  to  simplify  the  first  two  equations.  Making  the  substitutions 


,  _  1  ,  _  1  ,  1 
(l i  ?  £*2  t  2  5  ^3 


l  +  xz 


i+r 


l  +  z 


2  ’ 


(5.18) 


and  multiplying  the  first  two  equations  in  (5 . 1 7)  by  (1  +  x2  )(1  +  y 2  )(1  +  z2 )  simplifies 

these  equations  so  that  they  no  longer  depend  on  the  square  root  of  the  variables  we  are 

trying  to  solve  for.  The  three  equations  become 

((a2  -1  )y2  +  2abcos(0a  -0b)yz  +  b2z2  -c2)x2 
+(2bccos(0b  -  0c)(z 2  -  \)y  +  2accos(#fl  - 0c)(y2  -l)z)x 
-2 ab  cos {6a  -  0b )yz  +  (- b 2  +  c2z2  )y2  +  (1  -a2)z2  =  0 
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((/2  -  l)z2  +  2 df  cos (0d  -  0f  )yz  +  d2z2-g2  )x2 

+(2 dg  cos(0d  -  0g  )(y2  - 1  )z  +  2fg  cos(0f  -  0g  )(z2  -  l)y)x 

-2  df  cos  (0d -6>f)yz  +  (l-f2+g2z2)y2-d2z2=0  (5.19) 

1  1  1 

- 7  - 7  - 7  ~  P 

l  +  x  1  +  y  1  +  z 

Unfortunately,  the  first  two  equations  are  quadratic  in  x,  y,  and  z,  and  it  is  not  yet  known 
if  it  is  possible  to  find  an  analytic  solution  to  these  three  equations,  even  with  the 
additional  constraints  that  come  from  replacing  a-g  with  the  most  general  values  from  a 
SU(3)  matrix.  Additional  research  in  this  area  is  left  for  future  work,  and  I  have  focused 
on  a  specific  SU(3)  matrix  for  the  remainder  of  this  thesis. 


6  Diffusion  Equation 


6.1  Analytic  treatment 

The  diffusion  equation  can  be  modeled  if  the  SU(3)  matrices  inside  the  collision 
function  are  equal  to 


-in;  16 


U(3)  =  ■ 


V3 


( eiinn  x  x 

1  eillzn  1 


ilnl'i 


(6.1) 


Inserting  this  into  (5.12)  and  (5.13)  gives  the  collision  functions 


^i=| [-2 A  +  P2  +  A  - 2(2 a  -  1)VA(1-A)A(1-A)  + 

(2 A  -  1)VaO-A)a(1-A)  +  (2 A  - iWa(1-A)a(1“A) ] 


q2  =  }La  - 2 A  +  A  +  (2 A  - lW^(1-A)A(1-A)  + 

-2(2 a  - 1) V A  (!  -  A ) A  0  _  A )  +  (2 A  -  WaC1"  a)a(1_  A) ] 
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Q3  =  T [ P\  +  Pi  - 2 Pi  +  (2 Pi  - 0 (1  -  V2 )p3 0 -  Ji )  + 

i -  i -  (6-2) 

(2p2  - 1)7  a  (1  -  a  )  p3  (1- As)  -  2(2ft  -  !)VA  0  “  P\  )Pi  (]~Pi)] 
where,  as  expected,  Qj  +  Q2  +  Q3  =  0 . 

Running  the  numerical  simulation  presented  in  section  6.2  suggests  that  the 
equilibrium  values  dm ,  are  d]  =  d2  =  d3  =  p  /  3  .  This  can  be  easily  verified  by  noting  that 
at  equilibrium,  the  collision  matrix  should  not  change  the  occupation  probabilities;  that  is 
C\Yeq)  =  \Yeq)-  ThuS 


f  V3e*W6  0 


e~i7l/6 


0  1 
0  1 


0  0  0 

1  1  0 

2W3  j  o 

1  eiln‘l  0 


0  0  0 
0  0  0 
0  0  0 
0  0  0 


■sJdld2d3 
yjdld2{\  —  dz) 
y]dl  (1  -  d2  )d3 
yj{\-d{)d2d3 


yjdxd2d3 
7^2  (I- ^3  ) 
■yjdl(l  —  d2)d3 
y/O-  —  d{)d2d3 


(6.3) 


V3  0 


0 

0 


V 


0 


0  0 
0  0 
0  0 
0  0 


0  en’lz  1  1 

0  1  en’n  1 

0  1  1  e,2’n 

0  0  0  0 


0 

0 

0 

Set’16, 


vV(i-A)(i-^2)(i-^3)y 


V(i-</,)(i-rf2K 

vV(i-A)(i-^2)(i-^3), 


If  we  insert  dx  =  d2  =  d3  =  pi  3  into  the  above  equation,  the  expression  is  reduced 
to  the  following  equations 


V|2/3  f/^2/3 

V  3  J  v  3  y 

-ini  6 

p-  (2  +  e'>2/3 )  (d 

73  3' 


1/  \ 

1/  \ 

1  1-  — 

=  -J 

li-£ 

ll  3  J 

1  3  y 

3  y 

-(2  +  e'"2/3)  1— £ 

V  3  y 

\  2/3 


73 


V  3 


\-P 


\  2/3 


l-£ 


(6.4) 


which  of  course  can  be  further  simplified  to  the  identities  1  =  1  and  -£L-^—  (2  +  e  )  =  1 . 
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This  verifies  that  the  diffusion  equation  collision  matrix  has  no  effect  on  the  local  state 
when  d{  =  d2  =  d3  =  p  /  3  .  Thus,  this  must  be  the  equilibrium  condition. 

After  finding  the  equilibrium  condition,  the  next  step  is  to  perform  the  Chapman- 
Enskog  expansion  around  local  equilibrium.  From  section  4.3.2,  we  know  that  the 
equations 


2 

s  t — —+eem 


dt 


■*-.+s>ej°s£L+SeU‘!r*L+ (V)=n. 


dx 


dx 


2  dx 


(6.5) 


and 


ele 


— \  =  ej\5 /(0) 
dx  1 


(6.6) 


are  the  Taylor  series  expansion  and  first  moment  equation  of  the  FQLBE  respectively. 

The  vectors  in  (6.6)  are  now 

' g y(oO 

and  |  ,?/«")  =  Sf,m  (6.7) 

and  the  matrices  are 


8dm\_ 

(  ddl  \ 
dx 

del  2 

dx  / 

dx 

dd3 

J 

e  = 


(  ea{ 

anj 

Sfij  \ 

(  1  1  1  \ 

1  2  2 

f\ 

0 

0" 

df 

df2 

df 

0 

0 

0 

and 

J  = 

dCl2 

df 

dn2 

df2 

dQ2 

df3 

= 

1-11 

2  1  2 

1° 

0 

-b 

dn3 

dn3 

dQ  3 

1  1  -l 

V  2  2  V 

V  df 

dfi 

df  J 

fl,2,3=P/3 

(6.8) 


are 


Once  again  J  is  singular.  The  eigenvalues  and  left  and  right  eigenvectors  of  J 


\=- 


-1  2) 


E!  = 


0 

v1; 
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h=-- 

2 

(«2|  =  |(-1  2  "I)  |E2 

o 

II 

(^|  =  |(1  1  1)  |Ej 

r  \\ 


1 

v1/ 


(6.9) 


with  lengths  selected  so  that  (E.  \  Ej  j  =  Stj .  Since  ^^we  can  use  the  same  sort  of 
trick  Yepez  used  (described  in  section  4.3.2)  to  solve  for  It) :  multiply  both  sides  of 
the  first  moment  equation  (6.6)  by  J .  J  squared  is  equal  to 

W£,}(£,l  +  A|£2)(£2 1)2  =^(|£i)(£,|  +  |£2)(£2  |)(|£,)<£,  |  +  |£,)(£2|) 

=  ^(|£,)(£i|  +  |£2)(£2|)  (6.10) 

= v 

Thus  equation  (6.6)  can  be  solved  for  |t) /(0)^  as  follows. 

—\  =  j2\Sf(0)) 
etc/  1  ' 

f)=^l Sfm) 


Jte 

Jle 


te 


4 


a/' 

<9x  i 

ah 

dx  I 


implies 

=  T1|^/(0)) 

=  |j/(0)) 


(6.11) 


Inserting  this  solution  into  (6.5)  along  with  the  equilibrium  condition  dm=  p/3 
produces 

\  i  i 

-  +  0(e3)  =  Qm .  (6.12) 


2  1  ddm  2  2  „2  1  ddn 

£  T - -  +  £  etc  m 


3  dt 


3  etc  18 


dx 
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Summing  these  equation  over  m  with  e1=\,  e2  =  0,  and  e3  =  -1  produces  the 


diffusion  equation 


dp  _  r  ay 

dt  9  z  dx 2 


+  0(£3) 


(6.13) 


This  derivation  suggests  that  the  simulation  is  accurate  to  at  least  first  order  in 
time  and  second  order  in  space.  The  diffusion  coefficient  for  the  three  qubit  FQLGA  is 
equal  to  d  =  £2  /9z  where  £  is  the  lattice  spacing  and  z  is  the  time  step.  It  is  possible  to 
arbitrarily  change  either  £  or  z  to  adjust  the  diffusion  constant.  For  instance,  if  one 
wishes  the  diffusion  constant  to  be  1/9  m2/s,  then  £  is  one  meter  if  the  time  step  z  is 
defined  to  be  one  second. 


6.2  Numerical  treatment 

The  results  of  a  numerical  simulation  of  the  diffusion  equation  three  qubit 
FQLGA  model  are  presented  in  this  section.  All  numerical  simulations  were  run  in 
Mathematica  5.1  or  5.2  on  a  conventional  desktop  computer. 


6.2.1  Sum  of  Gaussian  and  sinusoid  initial  condition 

The  three  qubit  FQLGA  is  carried  out  using  circular  boundary  conditions — that  is, 
p(0,t )  =  p(L,t  )  where  L  is  the  length  of  the  lattice,  set  to  250/  for  this  simulation.  The 
initial  condition  was 


1  .  ( 2 7ZX^ 

p(x,  0)  =  —  sin 
6 


v  ^  J 


+lM)  +ii 

10  60 


(6.14) 


The  solution  to  the  diffusion  equation  with  circular  boundary  conditions  is 
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(6.15) 


(6.16) 


The  derivation  of  this  solution  is  presented  in  Appendix  B.  From  this  solution, 
one  should  expect  the  simulation  to  exhibit  a  roughly  exponential  decay  from  the  initial 
condition  to  the  average  of  the  initial  state  D0 ,  where  the  exponential  decay  constant  of 


higher  frequency  terms  (in  space)  is  larger  than  that  of  low  frequency  terms.  The  qubits 
are  each  initialized  to  be  p(z,  0)  /  3  so  that  they  will  be  close  to  local  equilibrium  after  the 


first  and  subsequent  time  steps.  The  results  of  the  simulation  are  presented  in  Figure  13. 

The  FQLGA  results  are  shown  in  gray  while  the  black  line  is  the  sum  of  the  first 
14  terms  of  the  analytic  solution  (6.15).  The  sum  is  cut  off  after  14  terms  because  the 
coefficients  Cm  and  Dm  for  m  >  1 4  are  less  than  the  negligible  value  10'10.  The  average 


percent  error  between  the  simulation  and  the  analytic  solution  is  shown  in  Figure  14  and 
Figure  16.  The  maximum  percent  errors  are  shown  in  Figure  15  and  Figure  17.  From 
these  figures,  one  can  see  that  the  magnitude  of  these  errors  oscillates  within  the  first  ten 
time  steps,  after  which  they  begin  to  steadily  decrease  as  the  simulation  approaches 
equilibrium. 
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Figure  13.  The  results  of  the  FQLGA  simulation  with  circular  boundary  conditions  for  a  number  of  time 
steps.  The  simulation  lattice  points  make  up  the  thick  gray  line  while  the  exact  solution  is  in  black.  The 
solution  decays  to  the  average  density  of  the  initial  condition. 
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Figure  14.  Plot  of  the  average  percent  errors  over  the  entire  lattice  for  the  first  20  time  steps.  Oscillations 
in  these  errors  are  apparent  at  the  beginning  of  the  simulation  but  disappear  after  about  ten  time  steps. 
After  this,  the  errors  decrease  as  shown  in  Figure  16. 
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Figure  15.  Plot  of  the  maximum  percent  errors  in  the  lattice  for  the  first  20  time  steps.  Oscillations  in  the 
maximum  percent  errors  are  also  apparent  at  the  beginning  of  the  simulation.  The  largest  percent  error 
during  the  simulation  occurs  after  the  first  time  step,  and  is  ~  0.14%. 
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Figure  16.  Plot  of  the  average  percent  errors  every  300  time  steps.  The  errors  decrease  as  the  simulation 
moves  closer  to  equilibrium. 


Figure  17.  Plot  of  the  maximum  percent  error  every  300  time  steps. 
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From  these  plots  we  see  that  the  maximum  percent  error  is  only  about  a  tenth  of  a 
percent,  while  the  largest  average  percent  error  is  .03  %,  indicating  that  the  simulation  is 
a  very  accurate  model  of  the  diffusion  equation  at  this  lattice  resolution.  Increasing  the 
resolution  (that  is,  increasing  the  number  of  lattice  sites)  only  improves  the  average 
percent  error,  as  is  clear  from  Figure  18.  This  is  a  log-log  plot  of  the  lattice  resolution 
verses  the  average  percent  error.  The  plot  was  made  by  running  identical  simulations  for 
15  time  steps  with  lattice  lengths  ranging  from  50 1  to  12800/ .  The  error  decreases  as 
percent  error ~ Sx2'01  where  Sx  =  £/L,  indicating  second  order  convergence  in  space. 


tf  L 

Figure  18.  Log-log  plot  of  the  average  percent  error  verses  the  lattice  resolution.  The  same  simulation  was 
run  with  lattice  lengths  ranging  from  50/  to  12800/ ,  each  evaluated  at  15  T  .  The  slope  of  the  best  fit 
solid  line  is  2.01,  indicating  second  order  convergence  in  space. 

6.2.2  Delta  function  initial  condition 

A  two  qubit  FQLGA  simulation  of  the  diffusion  equation  was  developed  by 
Yepez  [8].  This  algorithm  has  the  unfortunate  property  that  the  lattice  consists  of  two 
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interpenetrating  but  noninteracting  lattices.  This  is  due  to  the  fact  that  when  qubit  one 
from  lattice  site  m  streams  right,  qubit  two  of  lattice  site  m+ 1  streams  left.  Thus,  the  left 
streaming  qubit  reaches  lattice  site  m  without  ever  having  collided  with  the  right 
streaming  qubit,  which  is  now  on  lattice  site  m+ 1.  This  occurs  everywhere  in  the  lattice 
so  that  qubits  will  only  interact  with  those  from  every  other  lattice  site.  This  can  be 
demonstrated  if,  for  example,  the  simulation’s  initial  condition  is  a  delta  function1  as 


shown  in  the  left  column  of  Figure  19. 
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Figure  19.  Results  of  Yepez’s  two  qubit  diffusion  equation  simulation  with  a  delta  function  initial 
condition.  The  left  column  shows  the  results  of  the  unmodified  algorithm  with  qubits  one  and  two 
streaming  every  time  step.  The  sharp  spikes  in  these  pictures  are  due  to  the  noninteracting  nature  of  all  the 
qubits.  The  right  column  shows  a  modified  simulation  with  the  left  moving  qubits  streaming  every  even 
time  step  and  the  right  moving  qubits  streaming  every  odd  time  step.  The  modified  algorithm  corrects  the 
deficiency  in  the  algorithm.  This  figure  was  used  with  permission  from  [9]. 


1  This  is,  of  course,  a  Dirac  delta  function  with  height  equal  to  one. 
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There  are  two  methods  of  correcting  this  deficiency.  The  first,  taken  by  Yepez,  is 
to  stream  the  left  moving  qubits  every  even  time  step  and  the  right  moving  qubits  every 
odd  time  step.  The  results  of  this  process  are  shown  in  the  right  column  in  Figure  19.  An 
alternate  solution  to  this  problem  is  to  use  the  three  qubit  algorithm  I  have  developed. 

The  results  of  a  simulation  run  with  a  delta  function  initial  condition  are  presented 
in  Figure  20.  In  this  figure,  the  time  evolution  of  this  function  is  compared  to  the  analytic 
time  evolution  of  a  piecewise  defined  function.  This  piecewise  function  is  equal  to  zero 
everywhere  except  between  the  lattice  sites  adjacent  to  the  delta  function,  where  it  is 
triangular  in  shape  with  height  equal  to  one.  This  ensures  that  the  total  integrated  area  of 
the  piecewise  function  is  equal  to  the  total  density  of  the  lattice,  so  that  at  infinite  time, 
when  the  evolved  lattice  and  piecewise  function  have  reached  equilibrium,  the  constant 
lattice  density  will  equal  the  height  of  the  evolved  piecewise  function. 

The  three  qubit  algorithm  does  not  display  the  unusual  pattern  present  in  the 
unmodified  version  of  the  two  qubit  algorithm,  since  the  streaming  qubits  interact  via  the 
stationary  qubits.  However,  it  does  not  seem  that  this  method  is  superior  to  the  modified 
two  qubit  algorithm  for  two  reasons.  The  first  is  that  this  algorithm  can  be  much  more 
difficult  to  implement  experimentally  on  a  NMR  computer  because  increasing  the 
number  of  qubits  per  node  increases  the  complexity  of  the  system.  In  addition,  the 
diffusion  coefficient  for  the  three  qubit  algorithm  is  t 2  /  9t  ,  whereas  it  is  t 2  /  2 z  for  the 
two  qubit  algorithm.  This  means  that  for  equal  lattice  lengths  and  time  steps,  the  two 
qubit  algorithm  will  evolve  4.5  times  faster  than  the  three  qubit  algorithm. 

One  should  note  that  neither  diffusion  equation  algorithm  effectively  models  the 
evolution  of  a  large  gradient  function  with  poor  lattice  resolution,  such  as  a  delta 
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Figure  20.  Results  of  the  three  qubit  FQLGA  (black  data  points)  with  the  initial  condition  equal  to  zero 
everywhere  except  at  the  center  lattice  site  where  it  is  equal  to  one.  The  time  evolution  of  this  function  is 
compared  to  the  analytic  time  evolution  of  a  piecewise  defined  function  (solid  line),  with  an  integrated  area 
equal  to  the  total  density  of  the  FQLGA  simulation. 
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function.  This  is  because  the  “smallness  parameter”  £,  which  was  defined  to  be  on  the 
order  of  the  Knudsen  number,  is  no  longer  small.  The  Knudsen  number  is  equal  to  the 
mean  free  path  (the  lattice  spacing)  over  the  characteristic  length  (the  width  of  the  delta 
function).  Since  these  two  numbers  are  roughly  equal,  it  is  no  longer  appropriate  to 
neglect  the  terms  of  order  s 3  and  higher  in  the  equation 


dp_ 

dt 


e  ay 

9r  dx 2 


+  0(^3), 


(6.17) 


and  the  diffusion  equation  is  no  longer  a  good  approximation  of  the  governing  equation 
of  the  lattice.  It  is  not  surprising,  then,  that  the  agreement  between  the  algorithm  and  the 
analytic  solution  to  the  diffusion  equation  is  poor  at  the  beginning  of  the  simulation. 
Nevertheless,  at  times  after  about  16  r  the  Knudsen  number  becomes  smaller  and  the 
algorithm  begins  to  converge  to  the  analytic  solution  of  the  diffusion  equation. 

Another  case  where  the  FQLGA  can  be  expected  to  perform  poorly  when 
modeling  the  diffusion  equation  is  when  the  assumption 

/. =<r+«vr+^/:,+o(*3)  (6.i8) 


is  no  longer  valid — i.e.  if  8 /J;0)  /  dm  is  no  longer  on  the  order  of  s~  Kn.  In  this  case,  the 


Chapman-Enskog  expansion  will  fail  to  produce  an  equation  which  accurately  models  the 
lattice  in  the  continuum  limit.  For  example,  if  one  were  to  naively  choose  the  initial 
conditions  of  the  lattice  so  that  the  density  was  not  distributed  near  local  equilibrium  to 
begin  with,  then  8 /J;0)  /  dm  will  be  large  and  the  lattice  will  not  behave  in  a  manner  which 


mimics  the  diffusion  equation.  This  would  have  been  the  case,  for  instance,  if  the  density 
distribution  shown  in  Figure  13  had  initially  been  distributed  so  that  the  total  particle 
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density  at  every  lattice  site  was  located  in  a  single  qubit,  instead  of  being  spread  equally 
between  qubits  one,  two,  and  three.  Therefore,  to  model  fluid  dynamics  with  a  FQLGA  it 
is  essential  to  distribute  the  qubit  probabilities  near  local  equilibrium  when  a  lattice  is 
initialized. 

Returning  to  the  delta  function  simulation,  it  is  worth  pointing  out  that  this 
simulation  contains  the  most  extreme  gradient  possible  in  the  density  function  p .  For 
this  reason,  it  represents  a  strong  test  of  the  numerical  stability  of  both  the  two  and  three 
qubit  FQLGAs.  Both  algorithms  perform  well  under  these  conditions  due  to  the  unitary 
nature  and  structure  of  the  collision  operators.  Since  the  operators  are  unitary  and  block 
diagonal  to  preserve  the  probability  of  measuring  qubits  in  the  state  1 1) ,  the  total 

probability  of  finding  all  the  qubits  in  the  simulation  in  the  state  1 1)  is  preserved.  Thus, 

the  values  of  particle  density  will  remain  bound  and  not  “blow  up,”  and  are  in  this  sense 
numerically  stable. 

Another  common  definition  of  numerical  stability  concerns  the  accuracy  with 
which  an  algorithm  models  an  equation  [32],  It  is  clear  from  the  results  of  this  thesis  that 
the  three  qubit  Factorized  Quantum  Lattice  Gas  Algorithm  is,  in  this  sense,  a  numerically 
stable  model  of  the  diffusion  equation  when  the  Knudsen  is  much  less  than  one.  In 
addition,  the  results  presented  in  Figure  20  suggest  that  even  if  the  Knudsen  number  is 
initially  large,  the  algorithm  can  be  numerically  accurate  at  later  times  as  long  as  the 
initial  conditions  of  the  lattice  and  the  simulated  function  are  the  same,  and  the  integrated 
area  of  the  simulated  function  is  equal  to  the  average  density  of  the  FQLGA. 
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7  Conclusion 


This  thesis  was  written  in  support  of  AFRL’s  and  AFOSR’s  Quantum 
Computation  for  Physical  Modeling  basic  research  theme,  by  exploring  and  extending  a 
quantum  algorithm  designed  to  model  fluid  dynamics  using  a  practical  quantum 
computer.  To  date,  the  most  advanced  type  II  quantum  computer  prototype  uses  NMR 
technology.  Two  properties  of  the  algorithm  presented  in  this  paper  make  it  an  ideal  test 
case  for  the  modeling  utility  of  a  three  qubit  per  node  type  II  NMR  quantum  computer. 
The  first  is  that  there  is  a  maximum  of  three  qubits  that  the  algorithm  requires  to  be 
coherent  at  a  given  time.  This  allows  the  algorithm  to  be  run  on  a  set  of  three  qubit 
parallel  quantum  computers  connected  by  classical  communications  channels:  a  type  II 
quantum  computer.  The  second  property  that  makes  the  algorithm  ideal  for  a  NMR 
quantum  computer  is  that  information  is  stored  in  the  probability  coefficients  of  the  qubit 
basis  states.  Obtaining  these  probabilities  requires  an  ensemble  measurement  of  identical 
quantum  computers  running  the  same  algorithm,  which  is  precisely  the  sort  of 
measurement  used  in  a  NMR  machine.  Therefore,  this  algorithm  represents  a  good  test 
of  the  computational  capabilities  of  a  NMR  quantum  computer. 

This  thesis  extended  the  Factorized  Quantum  Lattice  Gas  Algorithm  from  two 
qubits  to  three,  in  an  effort  to  improve  the  algorithm  and  possibly  obtain  a  new 
macroscopic  governing  equation.  The  most  general  three  qubit  collision  operator  that 
preserves  particle  number  was  derived,  along  with  the  Quantum  Lattice  Boltzmann 
Equation  for  this  operator.  A  partial  derivation  of  the  governing  macroscopic  equation 
for  the  algorithm  in  the  continuum  limit  was  presented. 
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Difficulties  in  deriving  the  qubit  local  equilibrium  values  lead  me  to  consider  a 
more  specific  collision  operator.  For  this  operator,  the  governing  macroscopic  equation 
for  the  lattice  in  the  continuum  limit  was  the  diffusion  equation.  A  numerical  simulation 
of  this  algorithm  carried  out  on  a  conventional  desktop  computer  with  a  lattice  length 
L  =  25  Of  was  then  presented  and  compared  to  the  analytic  solution  of  the  one 
dimensional  diffusion  equation  with  circular  boundary  conditions.  The  simulation  and 
analytic  solution  matched  very  well — the  largest  average  percent  error  occurred  after  the 
first  time  step  and  was  0.03%.  Thereafter  the  error  decreased  as  the  simulation 
progressed.  Repeated  simulations  with  identical  initial  conditions  but  varying  lattice 
resolutions  revealed  that  the  simulation  possesses  second  order  convergence  in  space. 

Simulation  of  a  severely  under-resolved  gradient  (a  Dirac  delta  function)  was  also 
presented  to  check  for  numerical  instabilities.  No  numerical  overflows  occurred  and  the 
model  remained  stable  owing  to  the  collision  operator  being  unitary  and  block  diagonal, 
preserving  particle  number.  A  comparison  of  the  delta  function  evolution  using  the 
FQLGA,  to  the  analytic  time  evolution  of  a  piecewise  defined  function  revealed  that  the 
diffusion  equation  remained  an  accurate  governing  equation  for  the  FQLGA  in  the 
continuum  limit,  as  long  as  the  Knudsen  number  is  small.  In  addition,  it  was  observed 
that  the  Chapman-Enskog  expansion  will  not  result  in  a  valid  governing  PDE  unless  the 
ratio  of  the  lattice  deviation  from  local  equilibrium  to  the  local  equilibrium  density  is  on 
the  order  of  the  Knudsen  number. 
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Appendix  A.  SU(3)  Matrix 


A  unitary  matrix  may  be  generated  by  exponentiating  a  Hermitian  matrix. 
Exponentiating  a  Hermitian  quantum  mechanical  Hamiltonian  to  generate  a  unitary  time 
evolution  operator  is  a  well  known  example  of  this. 

The  collection  of  all  n  by  n  unitary  matrices  forms  a  closed  group  known  as  U(»). 
The  closed  subgroup  of  all  unitary  n  by  n  matrices  with  a  determinant  equal  to  one  is 
called  “special”  and  is  denoted  SU(«).  Obviously,  it  is  possible  to  create  members  of 
these  groups  by  exponentiating  Hermitian  matrices.  An  element  R  of  the  unitary  group  G 
may  be  written 

R  =  exp[/#S]  (A.l) 

where  the  Hermitian  S  is  called  a  generator  of  G  [31].  In  addition,  since  special  unitary 
matrices  have  a  determinant  equal  to  one,  the  generators  of  a  special  unitary  matrix  must 
be  traceless  [31]. 

det(R)  =  exp(tr(ln(R)))  =  exp  {id  tr(S))  =  1  =>  tr(S)  =  0  (A. 2) 


Since  SU(n)  is  a  closed  group,  multiplying  two  or  more  elements  within  the  group 
will  produce  another  element  of  that  group.  Multiplying  all  the  exponentiated  generators 
of  a  group  is  one  way  to  create  the  most  general  element  of  that  group.  For  instance,  a 
possible  choice  for  the  generators  of  SU(3)  are  the  Gell-Mann  matrices: 
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Notice  that  by  multiplying  each  of  these  matrices  by  a  free  parameter  and  summing  them 
together  one  will  obtain  the  most  general  three  by  three  traceless  Hermitian  matrix. 


The  exponentiated  Gell-Mann  matrices  are 
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Multiplying  these  matrices  together  is  one  method  of  generating  the  most  general  SU(3) 
matrix.  This  matrix  is  very  complicated  and  not  worth  explicitly  writing  here. 

For  the  purposes  of  the  FQLGA,  this  matrix  can  be  simplified  a  bit  if  the  order  of 
multiplication  is 


SU(3)  =  R8R7R6R5R4R3R2R1 

=  R8R7M 


(A.5) 


where  M  is  the  product  of  R6  through  Rj .  Inserting  this  notation  into  equation  (5.4)  we 
see  that 
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since  the  matrices  Dmn,  R7 ,  and  R8  are  all  diagonal.  For  instance,  for  the  matrix  D] 
one  obtains 
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Therefore,  for  the  purposes  of  this  FQLGA  it  is  possible  to  replace  the  full  eight 
parameter  SU(3)  matrix  with  the  six  parameter  matrix  M. 

In  fact,  preliminary  investigation  suggests  that  it  may  be  possible  to  replace  M 
with  a  simpler  four  parameter  matrix,  which  we  shall  label  N.  This  can  best  be 
understood  by  considering  that  the  purpose  of  the  matrix  N  is  to  redistribute  probabilities 
among  three  basis  states.  As  an  example,  suppose  we  were  interested  in  swapping  the 
probabilities  between  the  first  and  second  states.  One  could  then  use  the  following  swap 
matrix  to  do  so 
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The  determinants  of  all  the  root  matrices  are  complex  numbers  with  magnitude  one, 
meaning  they  are  all  unitary  as  opposed  to  special  unitary.  They  can  easily  be  made  a 
member  of  SU(3)  by  multiplying  each  matrix  by  its  determinate  to  the  -1  power. 
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However,  this  step  complicates  the  matrices  and  has  no  effect  on  the  final  results  of  the 
algorithm,  so  it  is  unnecessary. 

Thus,  a  simpler  four  parameter  matrix  which  can  replace  the  eight  parameter 
SU(3)  may  be 


SwapfSwap^SwapjSwapJ 


■  e* (}  +  eia  - eip  +  e‘(a+P)  +  e,y (1  +  3e'“  +  eip  - ei{a+p) )) 
f +  {  <f  *  (i  -  eia  -  eip  -  ei(a+P)  +  e‘r  (1  ■ -  3eia  +  eip  +  ei(a+/i) )) 


(A.  11) 


f +  \  e -» (f  -  eia  +  ew  -  ei(a+P) )  f +  j  e~¥  (j  +  -  eip  +  ei(“+/?)  -  e'y  (1  +  3eia  +  eip  -  e‘{a+P) )) " 

f + 1  e*  (f  +  eia  +  ew  +  ei{a+P) )  f -  }  e*  (f  +  eto  +  ew  +  ei(“+/?)  +  e,y  (1  -  3eto  +  eip  +  e'<Q+/?> )) 

\  +  \e~i*{=±  +  eip+eir(\  +  eip)) 


where  the  matrix  is  written  on  two  lines  due  to  its  length. 
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Appendix  B.  Analytic  Solution  of  the  Diffusion  Equation 


The  diffusion  equation 


di  dx 2 


can  be  solved  assuming  a  separable  solution  p  =  X(x)T(t) .  Then  the  diffusion  equation 
becomes 

x-T-=X=-e  (b.2  • 

d  T  X 

after  dividing  (B.l)  by  dXT .  Note  that  since  the  left  hand  side  of  the  equation  depends 
only  on  t,  while  the  right  hand  side  depends  only  on  x,  both  sides  must  be  equal  to  a 
constant  to  be  true  for  all  t  and  x.  The  constant  is  labeled  -k2  in  anticipation  of  what 


follows. 


By  inspection,  the  solution  to  the  equation  T' IT  =  - dk 2  is 


T{t)  =  Ae~dk  ‘  +B 


where  A  and  B  are  constants.  The  solution  to  X"  /  X  =  - k 2  is 

X(x)  =  C  sin(Ax)  +  D  cos(Ax)  +  E  (B  .4) 

where  C,  D,  and  E  are  also  constants.  With  circular  boundary  conditions, 

X(0)  =  D  +  E  =  C  sin (kL)  +  D  cos  (kL)  +  E  =  X(L)  (B.5) 

which  can  only  be  true  if  k  =  2nm  /  L  for  integer  m.  Of  course,  the  most  general  solution 
for  X(x)  must  be  the  sum  of  (B.4)  over  all  possible  m,  so  that 


p  =  X(x)nt)  =  B  +  Yjed{ML)t  Cm  sin  +Dm  cos 


2  nmx 


2nmx 
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The  initial  condition  is  thus 


**.(»  =  J  +  Z  C,,si„P=  +A,cosp= 

m=l  \  \  L  J  L 


To  solve  for  the  constants,  we  use  identities 


L 

r 

(  2nnx\ 

COS  -  1 

J 

0 

V  L  J 

L 

r  . 

(  2nnx\ 

sm  - u 

j 

0 

V  L  J 

L 

r 

( 2nmx\ 

cos 

J 

0 

V  L  J 

L 

r  . 

( 2nmx\ 

sm  - 

j 

0 

V  L  J 

dx  =  —S„ 


dx  =  —Sm 
2 


Then  multiplying  both  sides  of  (B.7)  by  sin(2 nnx/  L )  ,  integrating  from  zero  to  L,  and 
using  identities  two  and  four  gives 


c„  =yjp(x,0)sii 


In  nx 
L 


Similar  steps  allow  one  to  solve  for  the  constants 


B  =  —  |  p(x,  0  )dx 


2 1  f 

Dm  =  —  \p(x,  0)  cos 

Li  v 


2  nmx 
L 


(B-10) 
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